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Abstract 

The  objectives  of  this  research  are  to  develop  analytical  and  computer-aided  design  tech¬ 
niques  for  monolithic  microwave  and  millimeter-wave  integrated  circuits  (MMIC  &  MIMIC)  and 
subsystems  and  to  design  and  fabricate  those  ICs.  Emphasis  was  placed  on  heterojunction-based 
devices,  especially  the  High  Electron  Mobility  Transistor  (1IEMT),  for  both  low-noise  and 
medium-power  microwave  and  millimeter-wave  applications.  Circuits  to  be  considered  include 
monolithic  low  noise  amplifiers,  power  amplifiers,  and  distributed  and  feeuoack  amplifiers. 
Interactive  computer  aided  design  programs  have  been  developed,  which  include  large-signal 
models  of  InP  MISFETs  and  InGaAs  HEMTs.  Further,  a  new  unconstrained  optimization 
algorithm-POSM  has  been  developed  and  implemented  in  the  general-purpose  Analysis  and 
Design  program  for  Integrated  Circuit  (ADIC)  for  assistance  in  the  design  of  large-signal  non¬ 
linear  circuits. 

We  have  developed  an  accurate  analytical  drain  current-voltage  characteristics  model  for 
HEMT  devices  based  on  the  nonlinear  charge-control  formulation  which  describes  the  variable 
offset  distance  of  the  Two-Dimensional  Electron  Gas  (2-DEG)  from  the  heterointerface  in  a  more 
consistent  manner  and  with  greater  accuracy.  The  simplicity  of  this  model  makes  it  well  suited 
for  Computer-Aided  Design  (CAD)  applications  in  the  analysis  and  design  of  nonlinear  HEMT 
devices  and  integrated  circuits.  In  addition,  we  have  developed  a  general-purpose  finite-element 
two-dimensional  semiconductor  device  simulation  program,  which  is  able  to  analyze  and  simu¬ 
late  various  device  structures  iv.iuriirg  homo-  and  hetero-junction  III- V  compound 
semiconductor-based  devices  with  arbiu  geometries. 

The  work  also  included  the  synthesis,  growth,  characterization  and  device  prototype 
development  employing  InxGai_xAs  ternary  compound  semiconductor  material.  Pseudomoiphic 
InGaAs  HEMTs  with  higher  currents  and  therefore,  improved  power  performance  in  the  high  fre¬ 
quency  range  are  suitable  for  millimeter  wave  integrated  circuit  applications.  Single  heterojunc¬ 
tion  In2Ga.gAs/GaAs  and  double  heterojunction  In  isGassAs/GaAs  HEMTs  of  1  pm  and  .25 
pm  gate  length  have  been  fabricated.  Currents  up  to  310  mA/mm  and  510  mA/mm  have  been 
measured  for  the  1  pm  single  and  double  heterojunction  devices.  The  0.25  pm  devices  had  a  fmax 
of  85  GHz.  Also  reported  are  a  new  design  technique  for  power  distributed  amplifiers  using 
large-signal  S-parameters  from  harmonic  balance  analysis  and  design  examples  of  distributed 
amplifiers  using  the  MESFET  and  HEMT  devices. 
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Due  to  their  superior  high  frequency  and  high  speed  capabilities,  III- V  compound 
semiconductor-based  devices  appear  to  be  very  attractive  for  microwave  and  millimeter  wave  cir¬ 
cuit  applications.  New  fabrication  techniques  such  as  Molecular  Beam  Epitaxy,  high-resolution 
Electron  Beam  Lithography,  ...  etc.  allow  novel  high  performance  devices  to  be  realized  in  the 
laboratory. 

Computer  aided  design  techniques  have  to  keep  abreast  of  the  fabrication  technology  and 
precedes  it,  if  possible,  because  the  whole  design  must  be  made  a  priori  and  the  circuit  perfor¬ 
mance  are  simulated  via  CAD  tools,  if  we  adopt  a  monolithic  integrated  circuit  approach.  The 
computer-aided  design  of  microwave  circuit  is  now  commonplace.  Even  though  several  commer¬ 
cial  CAD  packages  are  available,  they  are  originally  developed  for  the  analysis  of  linear  elec¬ 
tronic  circuit  in  the  low  frequency  region  and  with  limited  number  of  device  models  available. 
Also  those  models  are  valid  in  a  narrow  operation  range  and  are  therefore  limited  in  their  applica¬ 
bility  to  the  latest  GaAs  MMIC  and  MIMIC  technology. 

There  is  a  growing  demand  for  the  improvement  of  the  CAD  tools  available  for  the  MMIC 
design  and  their  further  extensions  to  the  applications  for  MIMIC  design.  In  order  to  analyze 
MMIC’s  up  to  the  millimeter-wave  region,  a  simulator  should  be  able  to  analyze  a  circuit  design 
in  terms  of  its  topological  and  physical  parameters  including  frequency  dependence.  The  built-in 
device  model  of  the  circuit  simulator  should  be  accurate  and  simple  in  expressions  even  the  non¬ 
linearity  of  the  device  characteristics  is  taken  into  account  so  that  the  circuit  simulator  could  ran 
efficiently. 

The  high  electron  mobility  transistor,  HEMT,  has  become  an  important  device  because  of 
it’s  high  gain  and  performance  at  high  frequencies.  The  pseudomorphic  HEMT,  based  on  the 
strained  layer  InxGaj_x  As  epitaxial  material,  in  particular,  has  potential  for  being  the  device  of 
choice  for  millimeter  wave  integrated  circuit  applications[l,2],  HEMTs  employing  InGaAs  have 
higher  frequency  performance  than  the  conventional  AlGaAs/GaAs  HEMT  because  of  the  higher 
electron  mobility  and  velocity  of  the  InGaAs  compared  with  GaAs.  The  research  presented  is 
concerned  with  the  growth  of  pseudomorphically  strained  InxGa(_xAs/GaAs  heterostracturcs  and 
the  tailoring  of  the  composition  and  layer  structure  for  HEMTs.  Power  performance  of  devices  at 
high  frequency  will  be  an  important  issue  for  MMIC’s.  Prototype  devices  of  1  pm  and  0.25  pm 
gate  length,  with  higher  currents  have  been  fabricated  and  evaluated,  showing  the  potential  of 
these  devices  to  meet  these  needs. 

The  distributed  amplifier  is  based  on  the  principles  of  the  artificial  transmission  lines  into 
which  the  input  and  output  capacitances  of  an  active  device  are  absorbed,  and  then  the  gain- 
bandwidth  product  of  the  distributed  amplifier  can  be  increased.  Using  the  pseudomorphic 
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InGaAs/GaAs  devices  reported  here,  several  distributed  amplifiers  with  very  broad  bandwidth 
were  designed  and  simulated  at  UCSD. 
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II.  Computer-Aided  Design 

With  the  increasing  complexity  in  modem  semiconductor  device  structures  and  integrated 
circuits,  the  use  of  computer  programs  have  proven  to  be  valuable  aids  in  the  design,  develop¬ 
ment  and  characterization  of  new  devices  and  integrated  circuits.  The  process  of  CAD  consists  of 
three  important  segments,  namely  : 

I  Componet.  modeling 

II  Circuit  analysis 

III  Optimization. 

We  will  present  the  research  results  on  these  subjects  in  the  following  sections. 

2.1  Device  Modeling 

Device  modeling  has  played  an  important  role  in  the  computer  analysis  and  simulation  pro¬ 
grams.  Device  models  can  be  considered  in  two  broad  categories:  analytical  and  numerical 
models.  At  present,  the  relative  importance  of  analytical  and  numerical  models  is  clearly  esta¬ 
blished:  the  former  are  especially  useful  in  circuit  simulation  programs,  where  computational 
efficiency  is  the  main  concern.  The  latter,  instead,  arc  being  increasingly  used  for  verification  of 
analytical  model  and  as  design  tools  to  predict  the  achievable  performance  of  a  new  device  which 
has  not  been  fabricated.  III-V  compound  semiconductor  devices  show  great  promise  for  use  in 
high-performance  integrated  circuits.  Accurate  analytical  modeling  and  numerical  simulation  of 
modem  devices  become  quite  involved  and  time  consuming.  Difficulties  in  the  modeling  of 
active  devices  have  limited  the  use  of  CAD  techniques  at  microwave  and  millimeter  wave  fre¬ 
quencies.  To  aid  in  understanding  the  physical  operation  of  these  devices  and  to  optimize  their 
design,  accurate  device  models,  especially  for  heterostructure-based  devices,  arc  required.  The 
aim  of  this  research  is  to  develop  simple  and  accurate  analytical  models  and  efficient  numerical 
models  for  III-V  compound  semiconductor  devices,  which  can  help  in  design,  characterization, 
and  optimization  of  high-frequency  microwave  and  high-speed  digital  integrated  circuits. 

2.1.1  Analytical  model 

A  CAD-oriented  analytical  model  usually  contains  several  parameters  which  can  be 
adjusted  by  fitting  experimental  data;  so  doing,  reasonable  accuracy  can  be  achieved.  However, 
the  inadequacies  of  the  model  may  be  compensated  by  data  fitting  process  and  those  parameters 
can  hardly  be  traced  back  to  their  physical  meaning.  Thus,  in  order  to  link  device  performance  to 
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the  process  parameters,  simple  models  based  on  device  physics  with  a  minimum  number  of  fitting 
parameters  are  more  desirable.  In  other  words,  a  device  model  suitable  for  MIMIC  design  and 
optimization  should  encompass  the  following  features  :  (1)  It  should  include  the  effects  of 
process-related  device  parameters  (dimensions,  material  parameters,  doping  profile,  channel 
thickness,  recess  depth,  etc.).  Device-circuit  interactions  play  an  important  role  and  can  not  be 
ignored.  (2)  Device  models  must  be  computationally  efficient  to  allow  for  simulation  of  circuits 
in  an  economic  fashion.  (3)  The  model  must  provide  the  circuit  designer  with  synthesis  capabil¬ 
ity.  (4)  The  device  model  must  have  nonlinear  capability.  From  this  point  of  view,  we  have 
improved  the  charge  control  models  for  HEMTs  and  MISFETs  by  taking  the  nonlinear  variation 
into  account,  and  which  resulted  in  more  accurate  drain  current-voltage  characteristics  models. 

A.  HEMT  (High  Electron  Mobility  Transistor) 

The  analytical  models  most  widely  used  for  characterizing  HEMT  performance  arc  based 
on  the  linear  charge  control  model  which  either  neglects  the  variation  of  Fermi  potential  with  the 
applied  bias  or  assumes  a  constant  correction  distance  to  account  for  the  quantization  effect  in  the 
direction  normal  to  the  heterointerface  plane.  We  have  developed  a  simple  and  more  accurate 
nonlinear  charge  control  model  [3]  (see  Appendix  A),  derived  from  the  triangular  potential  well 
approximation.  Based  on  the  analytical  nonlinear  charge  control  formulation,  we  further 
developed  an  analytical  HEMT  drain  current-voltage  characteristics  model  [4]  (see  Appendix  B>, 
with  which  the  gain  compression  phenomena  near  the  pinchoff  regime  and  the  high  parasitic 
MESFET  conduction  are  more  accurately  described.  Moreover,  the  simple  and  analytical  form  of 
the  model  expression  make  it  very  suitable  for  CAD  applications. 

B.  MISFEJ  ( Metal-Insulator-Semiconductor  Field-Effect  Transistor) 

As  compared  with  the  empirically  optimized  Si-SiC>2  interface,  compound  semiconductor 
I-S  interfaces  are  generally  characterized  by  high  density  of  interface  states  with  a  non-uniform 
distribution  within  the  energy  band  gap.  To  obtain  more  accurate  description  of  the  modulation 
of  surface  potential  by  external  bias,  a  simple  empirical  distribution  of  interface  states  within 
energy  band  gap  is  included  in  the  charge  control  formulation.  The  results  show  that  the  inclusion 
of  interface  states  distribution  profile  into  drain  I-V  characteristics  model  leading  to  a  more  accu¬ 
rate  description  of  output  characteristics  of  III-V  MISFETs  [5]  (see  Appendix  C). 

2.1.2  Numerical  model 

As  device  configuration  and  fabrication  process  steps  become  complicated,  the  analytical 
approach  is  limited  in  its  accuracy  for  describing  two-dimensional  structures,  such  as  recessed 
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gate  transistors,  or  small-scale  devices  where  fringing  fields  and  other  multi-dimensional  effects 
predominate.  Accordingly,  numerical  models  are  necessary,  and  which  can  serve  as  (1)  an  aid  for 
a  deeper  understanding  of  device  behavior,  (2)  an  aid  for  the  development  and  validation  of 
CAD-oriented  analytical  models,  (3)  a  predictive  design  tool  for  device  performance  optimiza¬ 
tion.  We  have  developed  a  simple  and  feasible  finite  element  discretization  method  [6]  (see 
Appendix  D)  consistent  with  the  exponentially  fitted  Scharfetter-Gummel  scheme  for  variation  of 
carrier  densities,  in  our  general-purpose  two-dimensional  Semiconductor  Device  Analysis  pro¬ 
gram  (SDA).  Recently,  we  have  extended  the  application  of  SDA  program  to  the  simulation  of 
III-V  compound  semiconductor  devices,  including  homojunction  and  heterojunction.  The  simula¬ 
tion  results  of  two  example  devices  of 

A.  HEMT  (High  Electron  Mobility  Transistor) 

B.  MESFET  ( Metal-Semiconductor  Field-Effect  Transistor) 

are  reported.  The  simulated  distributions  of  potential  and  carrier  confirmed  the  validity  of  the 
assumption  of  nonlinearity  in  the  charge  control  formulation  used  in  the  analytical  I-V  model. 

Two-Dimensional  Finite-Element  Models  for  HEMT  and  MESFET  Devices 

High  electron  mobility  transistors  are  becoming  increasingly  important  for  the  provision  of 
superior  high  gri..,  high  performance  devices  for  microwave  and  miUimeter-wavc  applications. 
To  increase  the  transconductance  gm  and  cutoff  frequency  flt  device  geometries  must  be  scaled. 
The  future  trend  is  to  reduce  the  minimum  feature  size  down  to  deep  submicrometer  dimensions. 
As  the  feature  size  is  scaled  to  this  range,  careful  consideration  of  short  channel  effects  and  hot 
carrier  effects  in  the  conduction  channel  is  crucial.  To  this  end.  a  better  understanding  and 
analysis  of  device  operation  is  necessary  for  the  successful  development  of  device  and  accurate 
prediction  of  its  performance.  Two-dimensional  numerical  device  models  arc  suitable  for  investi¬ 
gating  the  details  of  intrinsic  complex  nonlinear  characteristics  in  the  scaled  devices. 

To  date,  papers  published  on  full  two-dimensional  models  for  HEMT  devices  arc  very  few 
[7-11].  In  addition,  these  models  are  restricted  to  finite-difference  methods,  which  arc  not  appli¬ 
cable  for  simulating  irregular  structures  such  as  recessed  gate  structures  commonly  used  in  the 
enhancement-mode  HEMT  devices. 

Recently,  we  have  successfully  developed  a  general-purpose  semiconductor  device  simula¬ 
tor,  which  is  capable  of  simulating  various  device  structures  with  arbitrary  geometries,  based  on 
new  finite-element  discretization  employing  the  Scharfettcr-Gummcl  (S-G)  scheme  [6].  By 
means  of  this  device  simulator  we  have  developed  two-dimensional  numerical  models  for  HEMT 
and  MESFET  devices  following  the  formulation  for  hctcrostruclurc  in  [12], 
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To  demonstrate  the  effectiveness  and  accuracy  of  the  new  two-dimensional  model,  the 
simulation  results  of  a  1  p.m  AlGaAs/GaAs  depletion-mode  HEMT  device  are  first  reported.  The 
simulated  DC  current-voltage  characteristics  of  the  HEMT  device,  whose  geometric  structure  and 
simulation  mesh  are  shown  in  Fig.  1  respectively,  is  plotted  in  Fig.  2.  It  is  obvious  that  the  gate 
current  plays  an  important  role  in  the  forward  gate  bias  region,  where  in  the  low  drain  bias  region 
the  source  current  dominates  and  shows  a  thermionic  emission  current,  which  is  suppressed  in  the 
reverse  bias  region.  The  potential  and  electron  distribution  plots  are  shown  in  Fig.  3  and  Fig.  4. 
From  the  distribution  plots  we  can  see  the  existence  of  two-dimensional  electron  gas  at  the 
heterointerface  due  to  the  conduction  band  discontinuity,  and  its  being  pinchcd-off  by  the  drain 
bias  resulting  in  the  current  saturation.  The  next  simulated  device  is  an  MESFET  device  with  gate 
length  of  0.7  pm,  whose  device  dimensions  and  simulation  mesh  arc  shown  in  Fig.  5.  The  result¬ 
ing  DC  drain  I-V  characteristics  are  shown  in  Fig.  6.  Fig.  7  and  Fig.  8  illustrate  the  potential  and 
electron  density  in  the  device  at  various  drain  voltages,  respectively. 

In  these  new  models,  we  have  also  included  the  energy  equation  to  take  short-channel 
effects  and  hot  carrier  effects  into  account.  It  is  planned  to  extend  the  two-dimensional  numerical 
model  to  the  modeling  of  submicron  gate  length  InGaAs-based  pscudomorphic  HEMTs  and 
AIGaAs/lnGaAs 

2.2  Circuit  Simulation  and  Optimization 

As  to  softwares  for  microwave  and  millimeter-wave  circuit  design  applications,  we  have  an 
in-house  program  called  NODAL  for  the  analysis  and  optimization  of  linear  high-frequency  ana¬ 
log  circuits  and  the  commercial  package  TOUCHSTONE  from  EEsof  for  the  same  purposes.  To 
design  nonlinear  large-signal  circuits,  we  have  CADNON  program,  which  is  an  in-house 
general-purpose  time-domain  simulation  and  analysis  program  based  on  state-variable  method  for 
microwave  nonlinear  passive  and  active  circuits.  In  addition,  we  have  successfully  ported  the 
SPECTRE  frequency-domain  simulator  for  nonlinear  circuits  to  Apollo  workstations  in  our  CAD 
laboratory.  The  SPECTRE,  developed  by  U.  C.  Berkeley  CAD  group,  uses  the  harmonic  balance 
approach  to  find  the  large-signal  response  of  a  nonautonomous  nonlinear  circuit.  Since  we  have 
source  codes  of  these  simulators,  we  can  tailor  these  programs  for  our  uses  by  incorporating  our 
new  device  models. 

Circuit  optimization  is  an  important  stage  during  the  design.  Chen  [13]  has  developed  an 
effective  optimization  algorithm  suitable  for  integrated  circuit  design.  However,  the  resort  to 
linear  search  is  still  frequent  in  Chen’s  method  and  the  amount  of  iterations  is  huge.  We  have 
developed  a  new  unconstrained  optimization  algorithm  using  the  Pseudo  Objective  Function  Sub¬ 
stitution  method  (POSM)  [I4|  (sec  Appendix  E).  This  algorithm,  requiring  neither  derivative 
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calculation  nor  linear  search  step,  substitutes  the  objective  function  by  a  second-order  approxi¬ 
mated  formulation,  which  enhances  the  convergence  rate  substantially.  The  algorithm  has  been 
implemented  in  the  general-purpose  analysis  and  design  program  (ADIC-2.C)  for  integrated  cir¬ 
cuit. 
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III.  Device  Fabrication 

3.1  Material  and  Device  Process 

In  addition  to  the  higher  electron  mobility  and  velocity  of  InGaAs,  the  advantages  for  using 
L.e  InGaAs  pseudomorphic  HEMT  are  that  the  use  of  AlGaAs  can  be  avoided  or  at  least  used 
with  a  low  A1  percentage  (less  than  20%).  Because  InGaAs  has  a  smaller  bandgap  than  GaAs,  the 
AEc,  conduction  band  edge  discontinuity,  is  comparable  or  greater  than  that  for  AlGaAs/GaAs 
HEMTs.  As  a  result,  the  problems  associated  with  trapping  centers  in  AlGaAs  of  higher  A1  per¬ 
centage,  needed  to  achieve  higher  AEc,  can  be  reduced.  Also,  the  wider  bandgap  of  GaAs  creates 
a  barrier  at  the  interface  between  the  InGaAs  and  GaAs  buffer  and  substrate.  This  serves  to 
improve  confinement  of  electrons  to  the  InGaAs,  reducing  parasitic  conduction  or  backdating 
through  the  buffer  or  substrate. 

In  order  to  achieve  higher  current  gain  cutoff  frequencies,  ft,  and  increased  currents,  which 
will  improve  power  performance  for  these  transistors,  it  is  desirable  to  increase  the  sheet  earner 
density,  ns,  without  degrading  the  electron  mobility  p.  This  can  be  achieved  through  increasing 
the  conduction  band  discontinuity,  AEC,  by  reducing  the  fundamental  bandgap  of  the 
InxGai_xAs,  with  larger  fractional  indium  concentrations,  x.  Increasing  AEC  also  decreases  the 
real  space  charge  transfer  back  into  the  lower  mobility  doped  layer.  An  increase  in  x  can  lead  to 
higher  electron  mobility,  but  this  may  be  offset  by  electron-electron  scattering  due  to  a  larger  ns. 
There  will  also  be  a  greater  lattice  mismatch,  resulting  in  a  larger  strain  so  the  thickness  must  be 
reduced  below  the  critical  thickness  for  the  generation  and  displacement  of  dislocations  at  the 
hetcrojunction. 

An  alternative  approach  to  increasing  the  sheet  carrier  concentration  is  to  employ  a  double 
hetcrojunction.  With  a  highly  doped  layer  both  above  and  below  the  active  channel,  charge 
transfer  can  occur  from  both  sides.  This  is  more  feasible  with  pseudomorphic  InGaAs  HEMTs 
because  the  quality  of  the  interface  between  InGaAs  on  top  of  GaAs  is  better  than  that  of  GaAs 
on  top  of  AlGaAs.  Results  will  be  presented  on  1  pm  and  0.25  pm  gate  length  single  heterojunc¬ 
tion  In2Ga8As/GaAs  HEMT  and  double  hcterojunctio.i  In isGa85As/GaAs  HEMT. 

The  InxGai_xAs  is  lattice  mismatched  to  GaAs,  so  its  thickness  and  composition  must  be 
chosen  so  that  it  is  elastically  distorted  relative  to  its  GaAs  substrate  thus  insuring  an  essentially 
defect  free  hetcrojunction  interface.  It  is  advantageous  that  HEMTs  be  developed  on  GaAs  sub¬ 
strates  because  they  are  of  significantly  better  quality  than  InP  substrates  and  arc  compatible  with 
current  MMIC  technology.  The  layer  thickness  of  the  In2Ga8As  is  only  135  angstroms.  HEMT 
structures  were  grown  with  25%  and  30%  In  fractions,  but  the  1pm  devices  that  were  fabricated 
had  inferior  performance  suggesting  that  on  a  GaAs  buffer,  using  InGaAs  with  In  percentage 
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much  above  20%  may  not  be  feasible  due  to  the  increased  lattice  mismatch  and  the  layer  thinness 
required  to  prevent  dislocations. 

The  electron  donor  layer  of  the  devices  is  n  type  GaAs,  avoiding  any  difficulties  that  may 
exist  with  using  AlGaAs.  The  doping  is  2.3xl018/cm3  from  Hall  effect  measurements.  The  layer 
structure  of  the  single  and  double  heterojunction  device  are  shown  in  Fig.  9(a)  and  9(b).  The  1 
pm  gate  length  devices  had  gate  widths  of  50  pm  and  150  pm.  The  device  fabrication  involves 
mesa  isolation,  AuGeNi  ohmic  contacts,  a  wet  chemical  gate  recess,  and  TiAu  gates.  The  0.25 
pm  gates  were  fabricated  on  wafer  by  direct  write  electron  beam  lithography  at  the  TRW  facili¬ 
ties  in  Redondo  Beach,  CA. 

3.2  Device  Results 

A  comparison  of  the  current  levels  of  the  1  pm  gate  single  heterojunction  In  oGa  gAs/GaAs 
device  and  the  double  heterojunction  In  15Ga85As/GaAs  device  with  the  single  heterojunction 
In  |5Ga  g5As/GaAs  device  done  in  earlier  work  are  favorable.  The  peak  current  of  'he 
In2Ga8As/GaAs  HEMT  was  310  mA/mm,  nearly  25%  better  than  the  250  mA/mm  found  for 
the  In  15Ga85As/GaAs  HEMT.  The  double  heterojunction  HEMT  had  peak  current  values  of  up 
to  510  mA/mm,  which  indicate  that  charge  transfer  is  occurring  across  both  interfaces.  These 
current  values  are  taken  near  the  peak  transconductance  values. 

To  achieve  higher  frequencies,  HEMTs  must  go  to  submicron  gate  lengths. 
In2Ga8As/GaAs  HEMT’s  of  0.25  pm  gate  length  have  length  have  been  fabricated  by  UCSD 
with  a  gm  of  up  to  500  mS/mm.  One  device  on  which  s-parameler  measurements  were  made  (up 
to  26  GHz)  had  an  extrapolated  maximum  frequency  of  oscillation,  fmax  of  85  GHz  and  current 
gain  cutoff  frequency,  ft,  of  40  GHz.  The  maximum  available  gain  (stability  factor  k  >  1  )  at  26 
GHz  was  10.6  dB.  An  example  of  the  well  behaved  dc  characteristics  of  these  devices  is  shown  in 
F'g.  10  and  a  plot  of  the  maximum  available  gain  is  presented  in  Fig.  11  One  feature  found  on 
many  of  the  devices  has  been  a  broad  transconductance  curve.  This  will  improve  the  operating 
range  of  the  devices.  Double  heterojunction  In.^Ga  85As/GaAs  HEMT’s  with  0.25  pm  gate 
have  exhibited  peak  currents  of  430  mA/mm,  and  a  gm  of  350  mS/mm.  A  fmax  of  75  GHz  and  ft 
of  50  GHz  were  extrapolated  from  measurements  on  one  device.  The  increase  in  f,  reflects  the 
higher  current  of  the  device.  The  dc  characteristics  are  shown  in  Fig.  12  and  the  plot  of  max¬ 
imum  available  gain  is  found  in  Fig.  13. 
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IV.  Circuit  Design 

4.1  Design  Approach  :  A  Systematic  Approach  to  Design  of  Distributed  Amplifier 


The  idea  behind  the  design  of  distributed  amplifier  is  to  avoid  the  input  and  output  capaci¬ 
tances  of  active  devices  being  frequency  limiting  parasitics  by  incorporating  these  capacitances 
into  the  per-unit-length  capacitance  of  the  artificial  transmission  lines  which  serially  connect  all 
the  input  ports  and  all  the  output  ports  together  in  a  distributed  manner.  The  topology  of  a  distri¬ 
buted  amplifier  is  shown  in  Fig.  14.  This  approach  allows  the  gain  of  active  device  to  be  paral¬ 
leled  without  the  corresponding  paralleling  of  the  device  input  and  output  capacitances,  and,  thus, 
enhances  the  overall  gain-bandwidth  product.  Previously  we  developed  analytical  and  computer 
aided  design  techniques  for  broadband  high  frequency  distributed  FET  amplifier  design,  and 
implemented  in  the  analysis  and  design  program  (DFETA)  [16].  The  approach  employed  in  the 
DFETA  program  can  be  summarized  as  follows.  For  simplicity,  a  unilateral  simplified  model  of 
FET  is  adopted.  The  gate  transmission  line  is  modeled  by  cascaded  lumped  low-pass  filter  sec¬ 
tions  which  have  the  same  characteristic  impedance  and  cutoff  frequency  as  that  of  the  gate 
transmission  line.  Because  of  the  very  small  drain  capacitance,  Cd,  of  FET,  the  corresponding 
characteristic  impedance  of  the  drain  transmission  line  is  too  large  for  a  realistic  design  of  output 
matching  circuit.  We  included  a  series  inductance  to  Cd  to  construct  the  equivalent  m-derived 
low-pass  filter  sections  which  emulated  the  drain  transmission  line  with  the  same  characteristic 
impedance  and  cutoff  frequency  of  the  drain  transmission  line.  The  m-derived  drain  transmission 
line  can  match  the  impedance  at  the  output  port  of  the  circuit  easily.  Gain  response  of  a  DFETA 
can  be  estimated  from  the  normalized  gain  curves  generated  by  the  DFETA  program  with  the 
normalization  parameters  (  xg  =  RgCgfc,  xd  =Rd  Cd  fc,  and  fn  =  f/fc)-  After  these  normalized 
design  parameters  are  adjusted  to  achieve  the  desired  or  optimal  flat  gain  response,  we  can  deter¬ 
mine  the  real  circuit  design  from  those  parameters. 

In  DFETA,  the  characteristic  impedances  of  gate  and  drain  transmission  line  and  the  DC 
voltage  gain  can  be  expressed  as  follows  : 
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fcg  :  the  cutoff  frequency  of  the  gate  transmission  line 
fd  :  the  cutoff  frequency  of  the  drain  transmission  line 
m  :  the  m-derived  parameter  (m  <  1) 

Cg  :  the  gate  capacitance  in  the  simplified  model 
Cd  :  the  drain  capacitance  in  the  simplified  model 
gm  :  the  transconductance  used  in  the  simplified  model 
N  :  the  number  of  FET  devices  used  in  the  circuit 
And  the  corresponding  lumped  element  values  can  be  expressed  as 
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The  lumped  inductance  can  be  introduced  by  using  microstrips  whose  size  can  be  determined  by 
the  following  formula  under  the  assumption  that  the  size  is  much  smaller  than  the  quarter  wave 
length  at  the  operating  frequency  [17], 

l  =  L-v\)pZo  (8) 

where 

1 :  the  microstrip  length 

L  :  the  value  of  inductance  Lg,  Ld,  or  Ldm 

Zq  :  the  characteristic  impedance  of  the  microstrip 

op  :  the  propagation  velocity  of  signal  in  the  microstrip 

Although  we  can  always  obtain  a  set  of  normalized  design  parameters  for  a  specified  or 
optimal  gain  response  from  the  simulation  results,  it  may  not  be  practical  for  those  components  to 
be  realized  in  a  monolithic  chip.  For  instance,  we  have  to  keep  the  active  devices  as  far  as  possi¬ 
ble  away  from  each  other,  e.g.  300  pm,  because  the  coupling  effects  among  those  active  devices 
is  a  big  concern  when  it  comes  to  the  monolithic  IC  design.  However,  from  equation  (l)-(3),  it  is 
apparent  that  increase  in  the  bandwidth  of  a  distributed  amplifier  will  cause  a  reduction  in  both 
the  characteristic  impedances  of  the  gate  and  drain  transmission  lines,  consequently  a  reduction 
in  the  length  of  microstrip.  In  order  to  avoid  unrealistic  design,  we  developed  a  systematic 
approach  to  designing  broadband  distributed  amplifier  which  takes  the  realizability  in  a 
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monolithic  chip  into  consideration. 

First,  we  derived  the  equation  governing  the  characteristic  impedance,  the  lumped  induc¬ 
tance  and  the  physical  microstrip  sizes  of  the  gate  and  drain  lines.  Some  example  cases  are 
derived  and  listed  in  Table  I  and  Table  II,  respectively,  where  the  1^,  values  are  omitted  because 
it  is  not  important  as  far  as  the  realizability  is  concerned. 

As  seen  from  the  data  listed  in  Table  I  and  Table  II,  there  are  restrictions  imposed  on  the 
design  of  distributed  amplifier  due  to  the  consideration  of  realizability.  The  results  indicate  that 
the  DC  gain  can  not  be  set  as  low  as  we  wish.  There  is  constraint  on  the  minimum  DC  gain.  So  is 
the  maximum  bandwidth  of  a  reasonable  design.  The  active  devices  employed  in  our  design  are 
0.15  pm  pseudomorphic  InGaAs/GaAs  HEMTs  which  were  fabricated  and  measured  by  TRW. 
Their  small-signal  equivalent  circuit  model  and  de-embedded  model  are  shown  in  Fig.  15(a)  and 
15(b),  respectively.  Based  on  the  de-embedded  model,  we  further  derive  its  simplified  unilateral 
model  as  shown  in  Fig.  15(c).  Referring  the  data  in  Table  I  and  II,  we  use  DFETA  simulation  pro¬ 
gram  to  evaluate  the  gain  performance  r'f  die  distributed  amplifier  design  where  we  used  four 
HEMTs.  The  best  results  of  the  distributed  amplifier  design  using  uniform  gate  and  drain 
transmission  lines  are  shown  in  Fig.  16.  We  further  incorporated  the  tapered  gate  and  drain  lines 
into  the  DFETA  simulation  program  in  order  to  enhance  design  flexibility  and  to  optimize  the 
gain  response,  in  which  we  set  the  source  and  load  impedances  of  the  distributed  amplifier  as  50 
ohms.  The  best  results  of  the  distributed  amplifier  with  tapered  lines  are  shown  in  Fig.  16.  To 
verify  the  design,  we  analyzed  the  distributed  amplifier  with  tapered  lines  by  the  TOUCHSTONE 
analysis/optimization  program  from  EEsof,  where  the  de-embedded  model,  instead  of  the 
simplified  model,  was  used.  The  results  are  shown  in  Fig.  17,  also  shown  are  the  results  after  par¬ 
tial  optimization.  The  results  show  the  gain  is  9.5  dB  and  the  bandwidth  is  40  GHz  of  the  lumped 
distributed  amplifier.  Finally,  we  converted  the  lumped  inductances  into  microstrip  design  on  a 
GaAs  substrate  with  thickness  of  100  pm.  The  analysis  results  of  the  initial  microstrip  design  are 
shown  in  Fig.  18,  also  shown  are  the  results  after  partial  optimization.  The  gain  is  about  8.5  dB 
and  the  bandwidth  is  42  GHz.  The  resultant  sizes  of  microstrips  and  other  element  values  are 
listed  in  Table  III. 

4.2  Power  Distributed  Amplifier  Design 

We  further  extended  the  application  of  the  systematic  approach  of  distributed  amplifier 
design  to  the  design  of  power  distributed  amplifier.  Because  the  operation  region  of  the  power 
device  may  include  the  nonlinear  region,  the  S-parameter  measured  under  the  large-signal  condi¬ 
tion  is  a  function  of  the  input  drive  power.  In  [18],  we  have  investigated  the  validity  of  applying 
conventional  large-signal  S-parameter  in  distributed  amplifier  design  by  harmonic  balance 
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analysis.  We  found  that  the  input  drive  power  level  plays  an  important  role  in  determining  the 
frequency  range  within  which  the  large-signal  S-parameter  can  be  applied  in  the  power  amplifier 
design  without  causing  severe  discrepancy. 

A  nonlinear  FET  model  was  used  in  our  initial  power  distributed  amplifier  design.  The 
model  is  shown  in  Fig.  19(a).  The  DC  characteristics  of  this  nonlinear  FET  model  are  shown  in 
Fig.  19(b).  The  DC  operating  point  of  the  FET  device  is  at  gate  bias  of  0.0  volt  and  drain  bias  4.0 
volts.  The  large-signal  S-parameter  beyond  6  GHz  derived  at  9  dBm  RF  input  power  drive  was 
shown  to  be  valid  in  [18].  The  corresponding  simplified  large-signal  model  is  shown  in  Fig.  20 
Following  the  systematic  approach  of  distributed  amplifier  design  proposed  in  the  preceding  sec¬ 
tion,  we  evaluated  the  power  performance  of  the  power  distributed  amplifier  by  DFETA  program 
using  the  simplified  large-signal  simplified  model.  The  results  are  shown  in  Fig.  21  with  a  com¬ 
parison  with  the  results  derived  by  Touchstone  using  the  large-signal  S-parameter  rather  than  the 
simplified  model.  The  complete  microstrip  design  of  the  power  distributed  amplifier  converted 
from  the  lumped  element  design,  including  the  bias  circuit,  is  shown  in  Fig.  22.  The  optimized 
output  performance  estimated  by  Touchstone  is  shown  in  Fig.  23.  To  evaluate  the  sensitivity  of 
the  power  performance  of  the  distributed  amplifier  design  to  the  input  RF  power  drive,  the  SPEC¬ 
TRE,  a  simulation  program  based  on  harmonic  balance  method,  was  used.  The  output  power 
responses  at  various  input  drive  power  levels  are  shown  in  Fig.  24,  where  the  curve  with  0  dBm 
input  power  drive  shows  the  power  gain  at  normal  or  small  signal  operation.  In  Fig.  24,  the 
curves  with  input  drive  power  lower  than  10  dBm  indicates  the  device  operates  in  the  linear 
region.  As  seen  from  Fig.  24,  there  is  a  significant  discrepancy  between  the  power  gains  below  6 
GHz  under  normal  operati  >  condition  estimated  by  SPECTRE  and  calculated  by  Touchstone.  It 
can  be  well  explained  that  below  6  GHz  the  large-signal  S-parameter  derived  at  9  dBm  input 
power  drive  are  no  longer  accurate  to  be  applied  in  the  broadband  amplifier  design,  since  their 
higher-order  harmonics  are  too  large  to  be  assumed  negligible  [18].  However,  higher  than  6  GHz 
the  power  gain  calculated  by  SPECTRE  is  almost  the  same  as  that  estimated  by  Touchstone.  The 
1  dB  compression  points  at  each  frequency  are  shown  in  Fig.  25,  which  indicates  the  maximally 
useful  large-signal  operation  of  a  power  distributed  amplifier.  These  preliminary  results 
confirmed  the  applicability  of  the  proposed  systematic  design  approach  of  distributed  amplifier  to 
the  power  distributed  amplifier  design  by  incorporating  the  method  of  deriving  the  accurate 
large-signal  S-paramctcr  by  harmonic  balance  analysis. 
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V.  Summary  and  Recommendations 

Analytical  and  numerical  device  models  for  III-V  compound  semiconductor  FETs  have 
been  developed.  By  consideration  of  the  nonlinearity  of  charge  control  mechanism,  the  analytical 
models  were  able  to  describe  the  general  shape  of  the  I-V  characteristics  more  accurately  without 
introducing  additional  fitting  parameters.  We  will  further  be  concerned  with  the  derivation  of  an 
appropriate  gate  capacitance-voltage  characteristics  model  taking  into  account  the  nonlinear  vari¬ 
ation  of  channel  charge  with  the  applied  electrical  field  and  its  frequency.  The  two-dimensional 
numerical  models  will  be  finished.  Then  the  short  channel  devices  will  be  investigated  to  further 
extend  the  applicability  of  the  numerical  model  to  simulation  of  the  hot  carrier  effects  which  can 
not  be  modeled  accurately  by  analytical  approach.  Also,  the  power  performance  of  those  devices 
will  be  evaluated. 

The  initial  results  are  very  encouraging  for  the  use  of  pseudomorphic  InGaAs  HEMT’s  in 
millimeter  wave  integrated  circuits.  With  In  ^Ga.gsAs/GaAs  HEMT  material  grown  at  UCSD, 
TRW  fabricated  0.15  pm  gate  length  HEMT  with  a  maximum  frequency  of  oscillation,  fmM,  of 
125  GHz  and  a  current  cutoff  frequency,  ft,  of  75  GHz.  In  2Ga  8As/GaAs  HEMT’s  of  0.25  pm 
gate  length  have  been  fabricated  by  UCSD  with  a  gm  of  500  mS/mm  and  fmax  of  up  to  88  GHz. 
Double  heterojunction  In  15Ga  85As/GaAs  HEMT’s  with  0.25  pm  gate  have  exhibited  a  current 
of  550  mA/mm,  a  gm  of  320  mS/mm  and  fmax  of  75  GHz.  The  next  step  will  be  fabrication  of 
wide  gate  devices  for  increasing  the  current  output  of  these  pseudimorphic  HEMT’s  for  use  as 
power  amplifiers.  Airbridges  will  be  used  for  the  interconnections. 

Work  in  design  method  for  distributed  amplifier  has  come  up  with  a  new  design  technique 
for  power  distributed  amplifiers  using  large-signal  S-parameter  derived  from  harmonic  balance 
analysis.  The  initial  designs  of  power  distributed  amplifiers  using  the  pseudomorphic  HEMT  dev¬ 
ices  have  shown  excellent  performance.  The  application  of  the  new  design  technique  to  the 
design  of  distributed  amplifier  operating  at  millimeter-wave  frequencies  is  under  way. 
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Abstract — A  simple  charge  control  model  of  the  two-dimensional 
electron  gas  (2-DEG)  of  HEMT’s,  which  explicitly  takes  into  account  the 
effective  distance  of  the  2-DEG  from  the  heterointerface,  has  been 
developed  for  nse  in  analytic  /-  V  and  C-  V  modeling.  In  this  model,  the 
Fermi  energy  level  versos  the  2-DEG  sheet  carrier  concentration  is 
represented  by  a  simplified  expression  derived  from  the  triangular 
potential  well  approximation  and  is  shown  to  be  dominated  by  terms  with 
different  functional  forms  in  two  distinct  operation  regions:  a  moderate 
carrier  concentration  region  and  a  subthreshold  region.  The  validity  of 
the  analytic  charge  control  model  is  supported  by  the  calculated  results  of 
self-consistent  quantum  mechanical  model. 

I.  Introduction 

HE  recent  development  of  the  HEMT,  employing  novel 
layer  structures  such  as  strained-layer  heterostructures, 
has  shown  a  promise  of  higher  transconductance  and  a  good 
power  capability  resulting  from  superior  electronic  transport 
properties  of  the  active  channel,  a  larger  sheet  carrier 
concentration,  and  better  confinement  at  the  heterointerface. 
These  enhanced  characteristics  lead  to  improved  device 
performance  which  demands  a  more  accurate  approximation 
of  the  charge  control  model  over  a  wider  operation  range  in 
order  to  obtain  accurate  HEMT  l-V  and  C-V  models  useful 
for  circuit  simulators.  The  analytic  models  most  widely  used 
for  characterizing  HEMT  performance  are  based  on  the  linear 
charge  control  model  [1J,  [2],  which  either  neglects  the 
variation  of  the  Fermi  potential  with  the  gate  bias,  or  for 
simplicity,  assumes  a  constant  correction  distance  in  the 
direction  normal  to  the  heterointerface  plane  accounting  for 
the  quantization  of  the  two-dimensional  electron  gas  (2-DEG) 
in  the  quasi-tri angular  potential  well.  These  assumptions  only 
hold  for  a  very  narrow  range  and  result  in  a  low  degree  of 
accuracy  in  the  model  when  device  behavior  is  characterized 
over  the  whole  operation  region. 

Recently,  Kola  et  al.  [3]  proposed  a  data  fitting  expression 
for  Fermi  level  versus  sheet  carrier  concentration  that  would 
take  into  consideration  the  modulation  effect  of  an  applied  bias 
on  the  Fermi  energy.  However,  the  expression  lacks  corres¬ 
ponding  physical  significance  and  requires  a  sophisticated  data 
fitting  process.  In  this  study,  we  developed  a  simple  yet  more 
accurate  Fermi  energy  expression,  derived  from  the  triangular 

Manuscript  received  July  7,  1988;  revised  September  26,  1988.  This  work 
w as  supported  by  AFOSR,  Boiling  Air  Force  Base,  DC,  under  Grant  AFOSR- 
86-0339  monitored  by  Dr.  G.  Witt. 

The  authors  are  with  the  Department  of  Electrical  and  Computer  Engineer¬ 
ing,  University  of  California,  San  Diego.  La  Jolla,  CA  92093. 

IEEE  Log  Number  8824829. 


potential  well  approximation,  which  yielded  a  more  rigorous 
analytic  charge  control  model  in  which  the  dependence  of  the 
effective  distance  of  the  2-DEG  from  the  heterojunction  on  the 
2-DEG  carrier  concentration  was  accurately  reflected. 

II.  Extended  Model 

To  obtain  an  exact  charge  control  formulation  of  the  2-DEG 
channel  in  HEMT  structures,  Poisson’s  equation  and  Schro- 
dinger’s  equation  would  need  to  be  self-consistently  solved. 
Unfortunately,  the  physical  calculations  are  too  involved  for 
use  in  analytic  device  modeling.  However,  an  approximation 
approach  based  on  the  linear  charge  control  concept  can  be 
applied  [4],  Additionally,  the  triangular  well  approximation, 
in  conjunction  with  Fermi  statistics,  was  found  to  give  a  fair 
description  of  the  2-DEG  carrier  concentration  dependence  on 
the  gate  voltage  if  an  appropriate  empirical  parameter  was 
incorporated  to  account  for  the  discrepancy  in  field  strength 
due  to  the  assumption  of  an  infinite  barrier  height  imposed  on 
the  junction  [5], 

Consider  a  HEMT  structure  with  an  n'-type  substrate, 
whose  energy  band  diagrams  are  shown  in  Fig.  1.  The  flat- 
band  voltage  is  the  gate  voltage  applied  simply  to  counterbal¬ 
ance  the  bending  of  the  quasi-Fermi  level  across  the  active 
channel  at  thermal  equilibrium  and  can  be  expressed  as 


where  <t>B  is  the  Schottky  barrier  height;  A Ec  is  the  conduction 
band  discontinuity;  V„  is  the  potential  difference  between  the 
Fermi  level  and  the  conduction  band  edge  in  the  neutral 
substrate  and  —  qND\{d\  -  de)2/ 2e(  is  the  potential  across 
the  fully  depleted  top  layer  under  flat-band  condition,  where 
C|,  du  and  NDl  are  the  permittivity,  thickness,  and  doping 
density  of  the  large-gap  semiconductor  layer,  respectively; 
and  d,  is  the  spacer  layer  thickness. 

When  the  gate  bias  Vc  is  higher  than  VFB,  thr  voltage  Vc  - 
VFB  can  be  partitioned  between  two  components;  V,  -  (V°b  - 
Vb),  the  voltage  drop  absorbed  by  the  top  depleted  layer,  and 
Vi,  the  band  bending  of  the  2-DEG  channel,  i.e., 

Vc-VFB=Vl  +  V2.  (2) 

Referring  to  the  energy  band  diagram,  the  band  bending  across 
the  2-DEG  channel  may  be  expressed  by 

V1=EF+Vn+Vc{x)  (3) 
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Fig.  1 .  Energy  band  diagram  of  a  HEMT  on  the  n  '  -type  substrate  in  the  2- 
DEG  charge  control  regime  by  a  Schottlcy  barrier,  (a)  Flat-band  condition, 
(b)  Thermal  equilibrium. 


where  EF  is  the  Fermi  potential  and  references  the  conduction 
band  edge  of  the  2-DEG  channel,  and  Vc(x)  is  the  channel 
potential  at  point  x  along  the  channel  defined  with  respect  to 
the  source  contact. 

According  to  the  triangular  potential  well  approximation 
and  Fermi-Dirac  distribution  the  relation  between  the  2-DEG 
carrier  concentration  n,  and  the  Fermi-level  can  be  established 
only  if  the  lower  and  the  first  excited  subbands  are  being 
considered  and  is  given  by 

£>=~  In  [exp  [(q/kT)y0  ’  nf3)] 

+  exp  [(q/kT)yx  •  /if3)]] 

[exp  [(q/kT)yx  ■  /if3)] 

-exp  [( q/kT)y0  •  /if3)]]2 
+  exp  [(qns/DkT)] 

•  exp  [(q/kT)( 70+71)  ’  «f3]^  j  (4) 

where  D  is  the  density  of  states  in  the  2-DEG,  and  y0  and  71 
are  material  parameters,  usually  estimated  from  experimental 
measurements  [1],  For  moderate  values  of  ns,  the  argument  of 
the  exponential  terms  is  small  and  a  series  expansion  of  these 
terms  converges.  For  qualitative  consideration,  we  end  the 
expansion  after  the  linear  term.  The  result  for  the  Fermi 
potential  takes  the  simplified  form 


( 


qn,  \ 
2DkT  ) 


D 

+  7  (7i~7o) 
4 


2 


•  <3  +  (7o  +  7.)  ‘  nsn-  (5) 

The  contribution  of  the  second  term  to  the  variation  of  Fermi 
energy  with  carrier  concentration  is  apparently  small  com¬ 


pared  to  other  terms.  Furthermore,  the  first  term  dominates 
only  in  the  subthre^hoftf  regtou.-In  the  normal  operation  region 
the  former  two  t^fms  can  be  considered^  be  a  constant  which 
can  be  evaluated  at  the  equilibrium  carrier  concentration 
Therefore,  we  can  approximate  £>  by  expression  of  the  form 

EF=EF0+y  •  /if3.  (6) 

By  charge  conservation  and  Gauss’  law,  further  manipula¬ 

tions  yield  the  2-DEG  charge  control  formulation  as 

Qn  =  qns=-AQc=-—^—[Vc-VT-Vc)  (7) 

di+Aa 

where  the  threshold  voltage  is  defined  as  VT  =  <j>B  -  (A Ec/q) 
+  Ep 0  —  V°b,  and  the  corresponding  effective  distance  of  the 
2-DEG  from  heterointerface  can  be  expressed  as 

Ad=—  •  n;u\  (8) 

Q 

This  analytic  expression  is  essentially  the  same  as  that  derived 
by  Stem  [6]  using  a  variational  method. 

HI.  Results  and  Discussions 

Fig.  2  shows  a  comparison  of  the  fitted  results  of  the  linear 
approximation  [2],  the  second-order  polynomial  expression 
[3],  and  (6)  with  an  exact  calculation  of  the  triangular  potential 
well  approximation  [1].  As  seen  in  Fig.  2,  (6)  has  a  better  fit 
with  the  exact  Ep  versus  ns  characteristics  than  other  approxi¬ 
mate  functions  in  all  regions  of  operation  of  interest,  where 
Epo  —  -  0.062  V  and  7  =  0.385  x  10*”.  The  correspond¬ 
ing  effective  distance  of  the  2-DEG  from  heterointerface  is 
also  compared  with  the  constant  effective  distance  used  in  [2], 
The  results  indicate  that  the  variation  of  the  effective  distance 
over  the  operation  range  significandy  affects  the  characteriza¬ 
tion  of  the  charge  control  of  the  2-DEG  in  HEMT  structure. 

To  further  confirm  the  effectiveness  and  consistency  of  this 
model,  the  numerical  calculations  of  electron  energy  levels  in 
a  GaAs/Gai_xAlxAs  heterostructure  by  exact  quantum  me¬ 
chanical  model  [7]  are  included  for  comparison  with  the 
results  of  our  model  and  are  illustrated  in  Fig.  3,  where  £>0  = 
-57.0  mV  and  7  =  0.298  x  10**.  The  slight  difference 
between  these  two  curves  in  the  low  carrier  concentration 
region  is  caused  mainly  by  the  subthreshold  effect  in  which  the 
variation  is  logarithmic  in  nature.  The  plot  of  corresponding 
effective  distances  shows  that  without  introducing  an  addi¬ 
tional  fitting  parameter  in  (8)  the  approximation  is  very  close 
to  the  numerical  results  of  the  quantum  mechanical  model  and 
indicates  that  the  extended  model  is  more  accurate  as  far  as  the 
nonlinearity  of  charge  control  of  the  2-DEG  is  concerned. 
Moloney  et  al.  [8]  have  also  given  a  similar  approximation  for 
the  Fermi  energy  versus  carrier  concentration  of  the  2-DEG 
while  modeling  the  C-V  characteristics  of  HEMT’s,  in  which 
a  functional  dependence  like  the  second  term  in  (5)  was 
assumed  dominant  at  room  temperature  rather  than  the  third 
one.  The  fits  based  on  this  formula,  as  shown  in  Fig.  3,  are 
fairly  good  at  first  sight  even  in  the  low  carrier  concentration 
region.  However,  the  corresponding  effective  distances  of  the 
2-DEG  from  the  heterojunction  deviate  considerably  from  the 
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Fig.  2.  Comparison  of  some  approximate  expressions  on  Fermi  energy 
versus  sheet  carrier  concentration  of  2-DEG  channel  at  300  K  and  the 
corresponding  effective  distance  of  2-DEG  from  the  heterointerface.  Note 
the  results  of  (6)  almost  overlap  those  calculated  from  the  triangular 
potential  well  approximation. 


S' 
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Fig.  3.  Comparison  of  the  exact  numerical  solution  and  the  approximate 
models  on  Fermi  potential  versus  2-DEG  density  in  a  GaAs/AlGaAs 
heterojunction  and  the  corresponding  effective  distance  versus  2-DEG 
concentration. 


results  of  the  exact  numerical  model.  This  comparison  reveals 
the  importance  of  the  extraction  of  a  proper  functional 
dependence  for  the  analytic  fitting  expression  from  a  detailed 
physical  analysis  to  ensure  the  consistency  of  the  formulation. 

In  conclusion,  a  better  general  agreement  with  the  exact 
calculation  of  the  quantum  mechanical  model  indicates  that  the 
extended  model  provides  a  realistic  description  of  the  2-DEG 
charge  control  of  a  HEMT  structure.  Further,  the  simplicity  of 
this  model  makes  it  well  suited  for  use  in  analytic  modeling  of 
HEMT’s.  It  has  been  successfully  applied  to  predicting  the 
I-  V  characteristics  of  pseudomorphic  InGaAs/GaAs  HEMT’s 
with  satisfactory  agreement  with  the  experimental  data  over  a 
wide  operation  regime. 
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Abstract 

Consideration  of  nonlinear  variation  of  the  sheet  carrier  concentration  of  the  Two- 
Dimensional  Electron  Gas  (2-DEG)  with  Fermi  potential  in  the  quasi-triangular  potential  quan¬ 
tum  well  of  the  HEMT  structure  has  led  to  a  bias-dependent  effective  offset  distance  of  the  2- 
DEG  from  the  heterointerface.  The  inclusion  of  the  variable  2-DEG  offset  distance  allows  model 
characterization  of  the  charge-control  mechanism  in  a  more  consistent  manner,  and  with  greater 
accuracy,  than  the  conventional  linear  charge -control  model  does.  Based  on  the  nonlinear 
charge-control  formulation,  we  developed  an  accurate  analytical  drain  current-voltage  charac¬ 
teristics  model  for  HEMT  device.  This  model  is  valid  over  a  very  wide  range  of  operation, 
extending  from  near-subthreshold  regime  to  high  parasitic  MESFET  conduction  regime.  This 
model  also  includes  the  broadening  effect  of  the  2-DEG  quantum  well  in  the  pinchoff  regime, 
providing  a  more  accurate  description  of  current  saturation  mechanism.  We  demonstrate  the 
effectiveness  and  accuracy  of  this  model  by  comparing  measured  and  modeled  dc  characteristics 
of  normally-on  as  well  as  normally-off  HEMT  devices.  Furthermore,  the  simple  analytical 
expressions  make  the  model  very  suitable  for  CAD  applications  in  the  analysis  and  design  of 
high-frequency  microwave  and  high-speed  digital  HEMT  devices  and  integrated  circuits. 

This  research  is  supported  by  AFOSR,  Boiling  Air  Force  Base,  DC,  under  Grant  No.  AFOSR-86-0339,  monitored  by  Dr. 

Gerald  Witt. 
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Introduction 

Since  it  was  demonstrated  in  1980,  the  High  Electron  Mobility  Transistor  (HEMT),  a  field- 
effect  transistor  that  takes  advantage  of  the  novel  properties  of  the  two-dimensional  electron  gas 
formed  at  the  heterojunction,  has  shown  very  promising  performance  for  botn  high-frequency 
microwave  and  high-speed  digital  circuit  applications.  In  order  to  optimize  the  performance  of 
HEMT  devices  and  circuits,  many  analytical  models  for  HEMT’s  based  on  a  simple  linear 
charge-control  model  have  been  proposed  [1-5].  These  models  are  similar  to  that  for  conventional 
MOSFET’s  with  only  the  gate  insulator  layer  being  replaced  by  a  fully  depleted  doped  semicon¬ 
ductor  layer.  These  models,  however,  show  a  limited  range  of  validity,  basically  because  of  the 
inaccurate  assumption  of  linear  charge-control  for  the  modulation  of  the  2-DEG  channel  charge 
by  the  external  bias.  Thus,  more  accurate  and  realistic  models  are  clearly  needed. 

The  inaccuracy  of  linear  charge-control  modeling  of  the  HEMT  operation  could  stem  from 
the  improper  estimation  of  the  effective  offset  distance  of  the  2-DEG  channel  from  the  heteroin¬ 
terface,  especially  in  the  pinchoff  regime.  Hughes  et  al.[ 6]  pointed  out  the  important  effects  of 
nonlinear  characteristics  of  the  charge-control  model  on  the  evaluation  of  transconductance.  Tak¬ 
ing  into  account  the  variable  2-DEG  offset  distance  from  the  heterojunction  by  numerical  calcula¬ 
tion,  they  obtained  a  fair  agreement  between  the  calculated  and  measured  transconductance 
characteristics  over  a  wide  operation  range.  However,  the  derivation  of  this  model  is  not  explicit 
in  terms  of  the  physical  parameters,  and  the  numerical  solution  is  too  complicated  for  applica¬ 
tions  in  the  analysis  and  design  of  integrated  circuits. 

Recently,  Kola  et  al.  [7]  proposed  a  data  fitting  expression  for  the  Fermi  potential  versus  2- 
DEG  sheet  carrier  concentration  that  would  take  into  consideration  the  variation  of  the  Fermi 
level  in  the  2-DEG  channel  with  the  externally  applied  bias.  In  order  to  avoid  a  sophisticated  data 
fitting  process  and  to  provide  a  direct  physical  insight  into  the  operation  of  the  device,  we  have 
developed  a  simple  but  more  accurate  nonlinear  charge-control  model[8],  which  is  derived 
directly  from  the  triangular  potential  well  approximation.  Based  on  the  analytical  nonlinear 
charge-control  model,  we  further  developed  an  analytical  HEMT  drain  current-voltage  charac¬ 
teristics  model,  in  which  the  gain  compression  phenomena  appeared  near  the  pinchoff  region  and 
in  the  high  parasitic  MESFET  conduction  region  are  more  accurately  described. 

Charge-Control  Model 

The  first  proposed  conventional  charge-control  model  [1]  derived  from  Poisson’s  equation 
for  HEMT’s  neglected  the  Fermi -potential  term.  The  resultant  model  (see  the  energy  band 
diagrams  in  Fig.  i )  shows  that  the  channel  carrier  concentration  depends  linearly  on  the  gate  bias. 
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where  ns  is  the  2-DEG  sheet  carrier  concentration;  £]  and  dj  are  the  permittivity  and  thickness  of 
the  electron  supply  layer  with  a  wider  energy-band  gap,  respectively;  q  is  the  electronic  charge; 
VG  is  the  gate-source  biased  voltage;  <j)B  is  the  barrier  height  of  the  Schottky-gate;  AEC  is  the 
conduction  band  discontinuity  at  heterojunction;  Npi  is  the  doping  concentration  in  the  electron 
supply  layer;  de  is  the  spacer  layer  thickness.  Note  that  subscript  1  refers  to  the  parameters  asso¬ 
ciated  with  the  electron  supply  layer  with  a  wider  energy-band  gap  and  subscript  2  for  the  2-DEG 
channel  layer. 

Drummond  et  al.[ 2]  extended  this  model  to  include  the  variation  of  the  2-DEG  Fermi 
potential  gas  with  the  sheet  carrier  concentration.  A  linear  approximation, 

VF  =  Vpo  +ans ,  (3) 


yields  the  linear  charge-control  formulation  used  by  most  of  the  existing  analytical  models, 

<*n’=TT77  tvG-VFo-VT],  (4) 

di  +Ad 

where  VFo  is  the  equilibrium  Fermi  potential;  and  Ad  is  the  corresponding  offset  distance  of  the 
2-DEG  from  heterointerface.  Ad  is  assumed  to  be  a  constant  value,  i.e., 

£ia  -  o 

Ad  =  —  =  80  A  (5) 

q 

This  model  was  further  modified  in  [6]  to  take  into  account  the  nonlinear  dependence  of  the 
average  position  of  the  2-DEG  channel,  Ad,  from  heterointerface  on  the  gate  bias.  The  2-DEG 
offset  distance  is  calculated  from  the  weighted  average  distance  of  all  energy  levels, 

Ad  = - - -  n  =  0,  1,2,  3,  •••  (6) 

ns 

where  n,n  and  dn  are  the  density  and  average  distance  from  the  heterointerface  of  carriers  in  the 
n-th  subband,  respectively.  The  total  carrier  concentration  ns  in  the  2-DEG  channel  is  a  sum  of  all 
occupancy  of  energy  states  at  various  energy  levels, 

n,  =  £nsn  n  =  0,  1, 2, 3,  •  •  •  (7) 


However,  the  computation  process  is  too  involved  to  be  useful  in  the  CAD  environment. 
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In  the  triangular  potential  well  approximation,  if  only  the  lower  and  the  first  excited  sub¬ 
bands  are  considered,  the  relation  between  the  2-DEG  sheet  concentration  ns  and  the  Fermi-level 
can  be  expressed  as 


Ep  =  (kT)  Ins 


,  kT 
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f  -j^rYi".5  ^Yon 
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+  eDkTekT 


(Yo+YiK5 
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where  D  is  the  density  of  states  of  the  2-DEG;  Yo  and  Yi  are  material  parameters  of  the  2-DEG 
channel  layer,  usually  estimated  from  results  of  resonance  experiments  [1].  Yo  and  Yi  link  the  2- 
DEG  carrier  density  with  the  longitudinal  quantized  energy  through  the  equations 

Eo  =  Yb  n?/3 

Ei=Yi  nP 


For  moderate  values  of  ns,  the  argument  of  the  exponential  terms  is  small,  and  a  series 
expansion  of  these  terms  in  (8)  converges.  We  can  then  approximate  the  Fermi  energy  Ep  by 

J_  _2 

Ep  =  (kT)  in(^~~)  +  ^~(Yi  ~  Yo)V  +  q(Yo  +  Yi )  ns3  (9) 


In  the  normal  operation  region  the  third  term  dominates  and  the  first  two  terms  can  be  con¬ 
sidered  to  be  constants  that  can  be  evaluated  at  the  equilibrium  carrier  concentration  ns0.  There¬ 
fore,  we  have  a  further  simplified  form 

2_ 

Ep=Epo+Yns3  (10) 

This  approximation  holds  as  long  as  the  device  is  not  operated  in  the  deep  subthreshold  region, 
where  the  quantization  effect  is  of  minor  importance  because  the  potential  well  broadens  remark¬ 
ably  and  the  subbands  are  closely  spaced.  An  advantage  of  using  (10)  is  that  only  two  adjustable 
parameters  are  involved.  They  can  be  determined  from  estimated  values  in  the  corresponding 
terms  in  (9).  This  requires  much  less  effort  when  compared  to  other  data  fitting  expressions.  Fig. 
2  shows  plots  of  the  fitted  results  of  the  linear  approximation^],  the  second-order  polynomial 
expression]?],  and  (10)  together  with  an  exact  solution  of  the  triangular  potential  well  approxima¬ 
tion  under  the  assumption  of  negligible  background  impurities  [1].  The  expression  in  (10)  clearly 
shows  a  better  description  of  Ep  versus  ns  characteristics  in  all  regions  of  interest,  where 
Ep0  =-0.062  eV  and  y  =  0.385  x  10“ 11  eV  m4/3. 

Substituting  (10)  into  (1)  and  rearranging  terms,  we  can  express  the  2-DEG  sheet  carrier 


concentration  as 
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where  V  is  the  channel  potential.  The  threshold  voltage  VTO,  the  linear  charge-control  coefficient 
p,  and  the  offset  distance  coefficient  £,  are  separately  defined  as 


AEC  cjNdi  •> 
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It  is  interesting  to  note  that,  from  comparing  (11)  and  (4),  the  corresponding  offset  distance 
of  the  2-DEG,  Ad,  from  the  heterojunction  is  proportional  to  n71/3,  which  is  essentially  the  same 
as  the  results  derived  by  Stem  [9]  using  a  variational  method.  In  fact,  the  key  feature  of  our  new 
nonlinear  charge-control  model  is  the  inclusion  of  the  nonlinear  dependence  of  the  offset  distance 
of  the  2-DEG  channel  on  the  2-DEG  sheet  carrier  concentration.  This  offset  distance  will  affect 
the  transconductance  performance  significantly[6].  Fig.  3  shows  plots  of  the  corresponding  offset 
distances  of  the  2-DEG  as  calculated  by  the  above  expression,  the  linear  charge-control  model, 
and  exact  quantum  mechanical  model  [10].  Without  additional  fitting  parameters  in  (11),  our 
approximation  of  the  effective  offset  distance  of  2-DEG  is  very  close  to  the  numerical  results  of 
the  quantum  mechanical  model,  with  a  maximum  discrepancy  of  less  than  10  A. 


Drain  I  -  V  Characteristics  Model 


/.  Empirical  Velocity  versus  Electrical  Field  Model 


Experimental  measurement  has  shown  the  velocity  saturation  phenomenon  in  III  -  V  com¬ 
pound  material  at  high  fields.  To  take  velocity  saturation  into  account,  the  following  functional 
form  is  commonly  employed  in  numerical  models  for  depicting  velocity  versus  electrical  field 
curve. 

pE  +  usal(-^  )" 

t)(E)  = - (15) 

1 +(#-)" 

Ec 

where  p.  is  the  low-field  mobility;  usai  is  the  saturation  velocity;  Ec  is  the  electrical  field 
corresponding  to  the  peak  velocity;  n  =  4  is  found  to  give  the  best  fit  to  the  experimental  data.  In 
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order  to  obtain  an  analytical  solution,  we  set  n  =  1  for  both  the  parasitic  MESFET  channel  and  the 
two-dimensional  electron  gas  channel,  i.e., 


v(E)  = 


iJsat 

<4  +  ^)E 


E<ES 


(16) 


=  \)sal  E  >  Es 


(17) 


where  the  electrical  field  Es  at  the  onset  of  velocity  saturation  is  determined  by  the  equation 

Es  =  -—  (18) 

2.  The  Threshold  of  the  Parasitic  MESFET  Conduction 


In  the  normal  operation  mode  of  HEMT’s,  the  depletion  of  the  doped  layer  is  from  the 
charge  transfer  and  the  Schottky  gate  depletion.  At  thermal  equilibrium  the  undepleted  channel 
width  of  the  parasitic  MESFET,  h,  as  shown  in  Fig.  1,  is 

2  (19) 


h  =  (d,  -  de)  - 


Not 


2et 

qNDi 
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where  Vn  is  the  potential  difference  between  the  conduction  band  edge  and  the  Fermi-energy  in 
the  neutral  carrier  supply  layer. 

The  threshold  voltage  Vc  of  the  parasitic  MESFET’s  conduction  can  be  easily  evaluated  by 
letting  h  =  0.  Then 


vc  =  4>B-vn- 


qNpj 

2£i 


(d!  -  de  - 


HsO  \2 
Not 


(20) 


3.  Intrinsic  HEMT  Operation  Mode  (VG  <  Vc) 

(a)  Linear  Region  (VD  <  Vsal) 

When  the  applied  drain  voltage  is  not  high  enough  to  accelerate  carriers  up  to  saturation 
velocity,  the  drain  current  can  be  evaluated  by  combining  (11)  and  (16)  and  using  a  one-step  sub¬ 
stitution  approximation  for  the  2-DEG  sheet  carrier  concentration.  We  then  obtain 
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ID  =  Wqn,u 

Wqpfe[VG-VTO-V] 
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Ec2 


where  W  is  the  gate  width;  L  is  the  nominal  gate  length;  t|  ( =  y  p2/3  )  is  the  nonlinear  charge- 
control  coefficient;  V  is  the  channel  potential;  Ec2  is  the  characteristic  electric  field  used  in  (16) 
for  the  2-DEG  case. 


Let  y  define  the  inverse  offset  distance  variable  as  a  function  of  channel  potential  for  the  2- 
DEG  channel  by 


y  =  [vG-vTO-v] 


1/3 


(22) 


The  boundary  conditions  on  the  2-DEG  channel  require  that  V(x)  equal  zero  at  x  =  0  and  VD  at  x 
=  L,  which  correspond  to 


yo  =  [VG-VTO]I/3  (23) 

yD  =  [VG  -  V-ro  -  VD]1/3  (24) 

Substituting  (22)  into  (21)  and  integrating  y  from  y0  to  yD  and  x  from  0  to  L,  we  obtain  the  drain 
2-DEG  current  below  saturation  given  by 
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where  A2  and  B2  are  defined  as 

Wp2qP 


A2  = 


L 

B2  =LEc2 


(26) 

(27) 


(b)  Saturation  Region  (VD  >  Vsat) 

When  the  2-DEG  carriers  are  accelerated  to  the  saturation  velocity  by  the  longitudinal 
electrical  field  parallel  to  the  current  flow,  the  2-DEG  channel  current  in  the  saturation  region  is 
expressed  in  the  equation 


Id  =  WqP 


(VG  -  Vto  -  Vsal) 


[1  +T|(Vg  -  VT0  -  Vsat)-I/3] 
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The  velocity  of  the  electrons  in  the  channel  reaches  saturation  at  the  point  x  =  Ls.  Further 
increase  in  drain  bias  will  increase  the  electrical  field  along  the  channel  and  cause  the  electrons  to 
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reach  velocity  saturation  at  a  point  closer  to  the  source  electrode,  i.e.,  reducing  the  effective 
length  Ls.  By  requiring  the  current  continuity  between  the  linear  regime  and  saturation  regime, 
we  arrive  at  the  expression  for  the  length  Ls  of  the  linear  channel  region 
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where  the  voltage  drop,  Vsal,  across  Ls  is 

vsaI  =  vG  —  Vjo  -y^ 


(30) 


4.  Current  Saturation  Mechanism 


To  determine  Ls  and  Vsat,  one  can  solve  the  two-dimensional  Poisson’s  equation  in  the 
pinchoft  region  by  using  self-consistent  approach  and  keep  the  lowest  space  harmonic  ol  the 
approximation  solution  along  the  boundary  [11],  The  approximation  yields  the  relation  between 
L,  and  Vsal  as 


2dsatEs2 

Vd  -  Vsat  =  - - sinh 
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n(L-Ls) 
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sat 


(31) 


Note  that  due  to  the  nonlinear  charge-control  mechanism,  the  effective  distance  of  the  saturated 
2-DEG  channel  from  the  gate  electrode  is  a  function  of  potential  at  Ls,  namely. 


dsat  d  1 


1  +T|  (  VG  -  VTO  ~  Vsat 


(32) 


As  illustrated  in  Fig.  4,  in  the  saturation  region  the  high  saturation  voltage  reduces  the  depth  of 
the  quasi-triangular  quantum  well,  and  the  2-DEG  channel  becomes  broader  at  me  drain  side. 
This  is  equivalent  to  an  increase  in  the  offset  distance  of  the  2-DEG  channel  from  the  heterojunc¬ 
tion.  This  phenomenon  could  be  observed  in  the  electron  density  distribution  plot  in  a  two- 
dimensional  numerical  HEMT  model,  e.g.,  [12,13]. 

The  value  of  ys  as  a  function  VG  and  VD  can  be  easily  found  by  combining  (29)-(32)  and 
using  a  simple  Newton’s  rooting  method  in  a  few  iterations  with  the  initial  value  set  as  yo-  Once 
y,  is  known,  the  drain  current  can  be  calculated  directly  from  (28). 


5.  The  Consideration  of  Parasitic  Components 


Once  the  gate  voltage  is  larger  than  Vc,  the  jitic  MESFET  channel  starts  conduction, 
i.e.  the  carrier  supplying  layer  is  not  fully  depleted,  leading  to  current  contribution  from  the  2- 
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DEG  channel  as  well  as  the  parasitic  MESFET  channel  underneath  the  Schottky  gate.  The  calcu¬ 
lation  of  parasitic  MESFET  conduction  current  is  very  similar  to  that  derived  in  [5],  The  detailed 
derivation  is  shown  in  the  Appendix. 

The  intrinsic  HEMT  model  does  not  explicitly  include  the  parasitic  resistances  Rs  and  Rd 
on  the  source  and  drain  sides.  To  take  these  parasitic  resistances  into  account,  we  can  calculate 
the  intrinsic  HEMT  bias  voltages  Vg  and  Vd  by  solving  the  equations 


Vd  =  Vds  -  ID  (Rs  +  Rd) 

(31) 

VG  =  Vgs-ID  Rs 

(32) 

where  Vds  and  Vgs  arc  the  terminal  voltages  at  the  drain  and  gate  electrodes  relative  to  the  source 
contact. 

Results  and  Discussion 

To  show  the  accuracy  and  effectiveness  of  our  nonlinear  charge-control  model,  we  compare 
our  calculated  drain  current-voltage  characteristics  with  experimental  measurements  of  a 
normally-off  HEMT  and  a  normally-on  HEMT  reported  by  Drummond  et  al.  [2]  and  Lee  et  al. 
[3],  respectively,  which  are  two  typical  samples  commonly  cited  by  the  linear  models  reported  in 
the  literature.  The  parameters  used  for  the  modeling  are  summarized  in  Table  1.  Fig.  5  shows  the 
I-V  curves  of  a  normally-off  HEMT  with  1  pm  gate  length  from  our  model  and  from  a  linear 
charge-control  model  in  conjunction  with  a  two-piece  velocity  model.  As  seen  in  "i-  5,  our 
modeled  results  are  in  good  agreement  with  the  experimental  data.  Here,  the  device  size,  Rs,  p, 
and  Ndi  are  measured  values.  Moreover,  we  do  not  need  to  add  a  resistance  parallel  to  the  chan¬ 
nel  as  in  [4],  in  order  to  fit  drain  I-V  curves  in  the  saturation  region  for  this  1  pm  gate  HEMT.  As 
indicated  in  the  electron  density  distribution  plots  of  a  two-dimensional  numerical  simulation 
[12,13],  the  carriers  in  the  saturation  region  (fewer  compared  to  those  in  the  linear  region)  are 
further  repelled  from  the  heterojunction  toward  the  substrate  by  the  gate  field.  When  fitting  the 
drain  current  curves  in  the  saturation  region,  because  of  the  basic  assumption  of  constant  offset 
distance,  the  linear  model  includes  only  the  current  contributed  by  the  carriers  within  the  constant 
effective  distance  from  heterointerface  and  leaves  out  others.  Therefore,  Rs  or  p  values  different 
from  measured  data  and  a  parallel  conductance  are  required  in  the  linear  model  to  take  into 
account  the  conduction  current  contributed  by  those  2-DEG  carriers  with  larger  offset  distance. 
Note  that  this  current  component  is  often  interpreted  as  the  substrate  current  or  leakage  current  in 
the  linear  model.  Since  the  channel  broadening  of  the  quantum  well  in  the  saturation  region  has 
been  included  in  the  nonlinear  model  by  using  the  bias  dependent  effective  offset  distance,  a 
more  accurate  estimation  of  carrier  concentration  in  the  2-DEG  channel  is  obtained.  This  results 
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in  a  more  accurate  prediction  of  the  drain  saturation  current  by  present  model  as  compared  to  that 
calculated  by  a  linear  model.  As  theoretical  estimation  or  experimental  measurement  of  satura¬ 
tion  velocity,  osat,  of  2-DEG  in  the  quantum  well  has  not  been  widely  reported,  we  use  the  simu¬ 
lation  values  calculated  by  Monte  Carlo  method  [14]. 

Fig.  6  shows  the  experimental  data  of  a  normally-on  HEMT,  our  calculation,  and  the 
modeled  results  by  a  linear  charge-control  model  using  a  three-piece  velocity  model.  The  parame¬ 
ter  values  used  for  modeling  are  also  listed  in  Table  1,  where  Rs  and  p  are  measured  values.  It 
can  be  seen  that  the  accuracy  of  the  derived  I-V  characteristics  is  quite  satisfactory  when  the  non¬ 
linear  modulation  effect  is  included  in  the  charge-control  formulation.  Fig.  7  shows  the  calculated 
and  measured  saturation  drain  current  as  a  function  of  gate  voltage.  In  the  near  subthreshold 
region,  our  results  agree  very  well  with  experimental  data,  due  to  the  more  accurate  nonlinear 
charge-control  description  embedded  in  the  drain  I-V  characteristics  model.  The  existing  linear 
models  usually  underestimate  the  2-DEG  sheet  carrier  concentration  and  fail  in  the  near  subthres¬ 
hold  region. 

Using  a  logarithm  functional  form  similar  to  that  of  the  first  term  in  (9),  which  proved  to  be 
dominant  in  the  subthreshold  regime,  a  refinement  in  modeling  the  nonlinear  charge-control 
could  further  extend  the  validity  of  the  present  model  into  the  subthreshold  region. 

Conclusion 

We  have  developed  an  analytical  HEMT  model  based  on  a  simple  nonlinear  charge-control 
formulation  that  explicitly  takes  into  account  the  variation  of  the  2-DEG  offset  distance  from 
heterointerface  as  a  function  of  bias.  The  validity  of  the  analytical  charge-control  model  is  sup¬ 
ported  by  the  results  of  a  self-consistent  quantum  mechanical  model.  The  I-V  characteristics  can 
be  e.  pressed  in  a  simple  analytical  form,  making  it  very  suitable  for  the  analysis  and  computer- 
aided  design  of  microwave  and  high-speed  HEMT  integrated  circuits.  In  addition,  although  the 
derived  model  is  still  in  an  approximate  form,  the  modeled  results  are  in  fair  agreement  with 
experimental  data.  The  model  can  therefore  be  applied  to  optimize  device  performance. 
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Appendix 

Inclusion  of  Parasitic  MESFET  Conduction  Model  (VG  ^  VG) 


Consider  the  case  where  the  MESFET  channel  is  pinched  off  at  x  =  Lm  by  the  gate  bias. 
The  channel  potential  at  this  point  is 

VM  =  Vp  -  (<t»B  -  Vn)  +  VG  ,  (Al) 

and  the  pinchoff  voltage  is  given  by 


Vr  = 


qNm  (d,  -  de  -  ~-)2 

^Dl 


2e< 


(A2) 


Within  the  conduction  channel,  wherever  the  channel  potential  is  lower  than  VM>  the  2- 
DEG  will  not  be  perturbed  by  the  gate  and  drain  voltages  and  will  stay  at  its  thermal  equilibrium 
values,  ns0,  which  can  be  estimated  self-consistentiy  by  the  equation 


El 

ns°  qdi  [1  +^n701/3] 


[4>B-Vn-Vp-VT0 


(A3) 


(a)  Linear  MESFET  Channel  (VD  <  VM) 


The  current  flowing  through  the  MESFET  channel  can  be  derived  as  in  the  conventional 
MESFET  model. 


Ii  = 


where  A,  and  Bi  are  defined  as 


,,  2  (4>b  -  Vn  -  VG  +  Vd)3/2  -  (<t>B  -  Vn  -  Vg)3/2 

VD  -  T  - 


(A4) 


Wqp^D]  ns0 

A]  = - ; - (di  -  de  -  -rj— ) 

L  Nm 

B]  =LECj 


(A5) 

(A6) 


In  the  case  of  VD  <  VM,  since  the  whole  2-DEG  channel  behaves  like  a  linear  resistor,  the 
current  flows  through  the  2-DEG  channel  is 

A2(4>B-Vn-Vp-VTO)VD 


h  = 


1  +  £  ns0  /3 


(A7) 
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The  total  drain  current  Id  is 
ID  =  Ii  +  I2 


(A8) 


(b)  Pinchoff  MESFET  Channel  (VM  <  VD  <  VMl) 

The  current  due  to  the  equilibrium  two-dimensional  electron  gas  can  be  calculated  as  in  the 
linear  region. 


I2  = 


LA2  (  <j>B  ~  Vn  ~  Vp  ~  VT0  )  Vm 
(LM  +  ^L)tl+^n701/3] 

e2 


(A9) 


The  current  flowing  through  the  MESFET  channel  can  be  evaluated  with  the  boundary  condition, 
V(Lm)  =  VM,  and  can  be  expressed  as 


MVM_|1Vp_<V-v.-vg>M] 

3  >/vT 


i.  = 


A2  [(pp  —  vn  -  Vp  —  v  J  I  I 

I2[l  +^n7o/3]  B2  Bi  J  m 


(A  10) 


Similar  to  the  derivation  in  the  normal  operation  mode,  the  2-DEG  current  can  be  explicitly 
expressed  as 


h 


-3A2  f  VMI4>B-Vn-Vp-VT0]  f  « 
n  +  3fl+^01/3] 

B2 


r 


+  ri6ln 


yp+o 

yM+ri 


(All) 


(c)  Saturated  MESFET  Channel  (Vsat  <  VD) 

In  this  operation  region,  formulation  of  the  current  flowing  through  the  MESFET  channel  is 
identical  to  (A10).  And  the  2-DEG  current  is  very  similar  to  that  in  the  pinchoff  region  and  can 
be  derived  as 


h  = 


-3A2L 


vM[<to-vn-vp-vn 


ro 


(Ls  + 


31 1  +  4nsd/3  i 


6 

-+z 

n=l 


-C2 


(-Dn 


i(6-n)(y"-ynM) 


+  rpln 


ys+~n 

yw+n 


(A  12) 
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Similarly,  the  variables  Ls  and  Vsat  can  be  solved  from  the  following  two  equations,  derived 
from  current  continuity  requirement  and  from  Poisson’s  equation  in  the  saturation  region  with 
fixed  boundary  d\ . 


(  Vq  -  Vxo  -  y,  )  Ec2  +  ES2  ys+rl 

—  —  -A  T,  X  A  ) 


EC2  EC2ES2 

vM[<5>B-vn-vp-vTO]  6,r  (-i)n 


1  3[1  +^n7o/3] 

2diEsl 

VD  -  Vsat  =  -  sinh 

K 


+  st 

n=l 


Jt(L-Ls) 


yi 


.(6-n)  (  n  _  n 


(y?  —  yo)I +ti°  in 


ys+n 


yo+n 


2d. 


|(A13) 
(A  14) 
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Figure  Captions 

Fig.  1  Energy  band  diagrams  of  a  normally-on  HEMT  on  an  n“-typc  substrate  (a)  at  thermal 
equilibrium  and  (b)  in  the  2-DEG  charge-control  regime. 

Fig.  2  Comparison  of  the  exact  numerical  solution  and  the  approximate  models  on  the  Fermi 
potential  versus  2-DEG  sheet  carrier  concentration  in  a  GaAs/AlGaAs  heterostructure  at 
300°  K.  Note  that  the  fitted  results  of  equation  (10)  are  almost  indistinguishable  from 
the  exact  calculations. 

Fig.  3  The  effective  offset  distance  of  the  2-DEG  from  the  heterointerface  as  derived  from  the 
results  shown  in  Fig.  2  by  both  the  linear  and  nonlinear  charge-control  models. 

Fig.  4  Schematic  diagram  shows  the  2-DEG  channel  spreads  out  toward  the  substrate  in  the 
saturation  regime  and  the  conduction  of  parasitic  MESFET.  The  decrease  of  the  2-DEG 
sheet  carrier  concentration  and  the  broadening  of  the  quantum  well  result  in  an  increase 
in  the  effective  offset  distance. 

Fig.  5  Measured  and  calculated  drain  current-voltage  characteristics  for  a  normally-off  HEMT 
[2].  Solid  line:  present  model.  Dots:  measured  points.  Dashed  line:  linear  charge-control 
model  with  a  two-piece  velocity  approximation. 

Fig.  6  Measured  and  modeled  drain  I-V  characteristics  of  a  normally-on  HEMT  [31.  Solid  line: 
present  model.  Dots:  measured  points.  Dashed  line:  linear  charge-control  model  with  a 
three-piece  velocity  approximation. 


Fig.  7  The  calculated  and  experimental  saturation  drain  current  of  a  normally-on  HEMT  [3].  as 
a  function  of  gate  bias  (VD  =  1.5  V).  Solid  line:  present  model.  Dots:  measured  points. 
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Table  1 

The  device  parameters  used  by  present  model  for  modeling 

the  HEMT’s 


DEG  Sheet  Carrier  Concentration  (m  2 ) 


quantum  mech.  model 


Fig. 


( mA ) 
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Analytical  and  Computer-Aided  Models  of 
InP-Based  MISFETs  and  Heterojunction  Devices  * 


Presented  at  Joint  NOSC/NRL  InP  Microwave/Millimeter  Wave  Technology  Workshop,  San  Diego,  CA,  January  25-26,  1989. 


Analytical  and  Computer-Aided  Models  of 
InP-Based  MISFETs  and  Heterojunction  Devices 


The  superior  properties  of  InP  material,  e.g.,  higher  peak  electron  drift  velocity,  thermal 
conductivity,  and  breakdown  field,  to  GaAs  have  made  it  an  alternative  for  high  performance 
applications  in  microwave  and  millimeter-wave  regimes  as  well  as  high-speed  digital  circuits. 
Recently,  a  high-efficiency  InP  MISFET  has  demonstrated  4.5  watts  output  power  with  4dB  gain 
and  46%  power-added  efficiency  at  9.7  GHz  by  Messick  et  al  [15].  These  impressive  results 
clearly  confinmed  the  promising  superior  performance  og  InP  MISFETs. 

The  main  concern  in  the  applications  of  I1I-V  MISFETs  has  been  the  reliability  of  output 
characteristics  of  the  devices,  which  is  mainly  attributed  to  the  variations  of  interfacial  properties 
of  the  gate  dielectric  layer  and  the  underlying  semiconductor  active  layer.  The  task  of  modeling 
output  characteristics  of  III- V  compound-based  MISFET  devices  has  become  complex  with  the 
possible  dominance  of  the  interfacial  properties  in  the  devices’  performance.  Much  more  atten¬ 
tion  should  be  paid  to  the  nonlinear  modulation  of  the  surface  potential  by  the  external  gate  vol¬ 
tage  due  to  the  presence  of  an  excessive  amount  of  interfacial  states,  since  the  accompanying  car¬ 
rier  trapping,  scattering,  and  recombination  could  have  altered  completely  the  charge  control  and 
transport  mechanisms,  and  consequently  the  device  output  characteristics. 

Because  of  these  factors,  we  have  developed  analytical  and  computer-aided  models  for 
depletion-mode  and  accumulation -mode  MISFETs  based  on  a  nonlinear  charge  control  model 
derived  from  semi -cm  pin  cal  surface  potential  formulation,  which  can  provide  us  not  only  an 
accurate  description  of  the  drain  I-V  characteristics,  but  also  a  comprehensive  study  of  the 
influence  of  interfacial  properties  on  the  output  performance  of  InP  MISFETs.  Key  aspects  of  the 
physics  of  this  device,  which  relate  to  charge  control,  carrier  trapping,  and  field-dependent  mobil¬ 
ity,  arc  modeled  in  this  study. 

In  order  to  further  verify  the  calculated  results  and  predict  device  performance  in  the  submi¬ 
cron  gate-length  regime,  we  have  developed  a  general-purpose  finite  dement  two-dimensional 
semiconductor  device  simulation  program,  whic'  is  able  to  analyze  and  simulate  various  device 
structures  including  homo-  and  hetero-junctions  III-V  compound  semiconductor-based  devices 
with  arbitrary  geometries.  Preliminary  simulation  results  of  a  1  pm  AlGaAs/GaAs  HEMT  device 
arc  reported.  It  is  planned  to  extend  the  two-dimensional  model  to  submicron  InP  MISFETs. 
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Abstract— The  finite  element  method  has  been  playing  an  important 
role  in  the  solution  of  many  engineering  problems.  To  make  the  appli¬ 
cation  of  the  finite  element  method  practical  in  semiconductor  device 
simulation,  we  have  applied  the  Scharfetter-Gummel  (S-G)  scheme  in 
conjunction  with  an  accurate  seven-point  Gaussian  Quadrature  rule  to 
the  assembly  of  the  finite-element  stiffness  matrices  and  right-hand  side 
vector  of  the  semiconductor  equations,  which  has  resulted  in  solutions 
with  high  accuracy  even  on  coarse  mesh  as  well  as  a  significant  speed¬ 
up  of  convergence  rate. 


I.  Introduction 

HE  Scharfetter-Gummel  scheme  [1]  has  been  widely 
used  in  the  discretization  of  the  continuity  equations 
and  has  shown  great  suitability  due  to  an  inherent  advan¬ 
tage  of  uniform  convergence.  However,  to  date  the  ap¬ 
plications  have  been  almost  restricted  to  the  finite  differ¬ 
ence  method  or  its  variants.  The  classical  finite  element 
method,  by  which  the  carrier  densities  are  usually  ap¬ 
proximated  in  low  order  polynomial  expansions  [2],  re¬ 
quires  an  excessively  fine  finite  element  mesh  in  order  to 
achieve  a  reasonable  approximation  error.  Additionally, 
it  suffers  a  prohibitive  convergence  rate  because  of  the 
large  integration  errors  incurred  due  to  the  inaccurate  low 
order  polynomial  approximations.  Therefore,  the  classi¬ 
cal  finite  element  method  is  inferior  to  its  finite  difference 
counterpart  employing  the  Scharfetter-Gummel  (S-G) 
scheme  in  which  the  functional  variation  of  carrier  den¬ 
sities  is  exponential  in  nature,  the  so-called  Box  Integra¬ 
tion  Method  (BIM-SG)  f3],  in  the  treatment  of  current 
continuity  equations. 
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Recently,  various  upwind  schemes  [4J  have  been  pro¬ 
posed  for  finite  difference  and  finite  element  to  achieve 
the  upwind  effect  in  order  to  avoid  the  spurious  oscillation 
in  solving  the  convection  dominated  flow  equations.  In 
essence,  the  optimal  upwind  formulation  in  one-dimen¬ 
sional  case  (without  source  terms)  can  reduce  to  a  form 
equivalent  to  the  S-G  scheme.  Nevertheless  the  flow  di¬ 
rection  not  always  being  along  the  grid  line,  it  is  difficult 
to  extend  these  formulations  to  two  or  three  dimensions 
in  a  straightforward  manner  [5]. 

A  simple  and  feasible  finite  element  method  consistent 
with  the  exponentially  fitted  scheme  has  not  been  devel¬ 
oped  for  two-dimensional  problems.  In  order  to  further 
exploit  the  promising  usefulness  of  the  S-G  scheme  in  the 
finite-element  method,  we  have  developed  in  our  general- 
purpose  two-dimensional  Semiconductor  Device  Analysis 
program  (SDA-1)  a  simple  and  accurate  finite  element 
discretization  method  (FEM-QSG).  where  the  S-G  scheme 
is  embedded  in  the  quadrature  of  finite  element  assembly 
of  Poisson  equation  and  continuity  equations.  The  new 
method  not  only  avoids  the  problem  of  flux  direction  in 
the  discretization  of  flow  term,  but  also  extends  the  ap¬ 
plication  of  the  S-G  scheme  to  the  discretization  of  source 
term  in  a  consistent  way. 

In  Section  II,  the  basic  governing  equations  in  SDA-1 
and  the  classical  finite  element  discretization  method  are 
reviewed.  Next,  a  new  finite  element  discretization 
method  employing  the  S-G  scheme  is  presented  in  Section 
III.  In  Section  IV,  some  typical  simulation  results  are  dis¬ 
cussed  and  are  compared  with  the  results  obtained  by  the 
classical  finite  element  method. 

II.  Formulation  of  Basic  Equations 

In  order  to  devise  a  general-purpose  device  simulator, 
which  is  capable  of  simulating  various  device  structures 
including  homojunctions  and  heterojunctions  with  arbi¬ 
trary  geometries  under  any  operation  condition,  we  based 
our  analysis  on  a  macroscopic  description  of  semiconduc¬ 
tors.  In  accordance  with  the  quasi-Fermi  level  concept  and 
in  terms  of  normalized  variables,  the  basic  semiconductor 
equations  governing  device's  operation  with  (yL  v,  w)  (  r 
=  e~°\  w  =  e *p)  as  dependent  variables  are  the  following 
[6].  (7]: 

V  •  (eY<£)  =  (n  -  p  -  No  +  N;)  ( 1 ) 
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V  •  (nae«  +  K)V»)  =  -  +  R  (2) 

V  •  (lxpe(~*  +  v’)V<j)  =  ^  +  R  (3) 

where  4>n  and  4>p  are  the  quasi-Fermi  potentials  for  elec¬ 
trons  and  holes;  v  and  w  are  the  Slotboom  variables;  e  is 
the  permittivity  of  semiconductor;  Np  and  N ~A  are  the 
concentrations  of  ionized  donor  and  acceptor;  pn  and  pp 
are  the  mobilities  of  electrons  and  holes;  R  is  the  net  re¬ 
combination-generation  rate  including  Shockley-Read- 
Hall  recombination,  Auger  recombination,  impact  ioni¬ 
zation  .  .  .  ;  n  and  p  are  the  electron  and  hole  densities 
and  are  given  by 

n  =  nire(*  +  v’~*")  (4) 


(  $  +  v*  —  ) 


p  =  nire 


($p  -  ^  +  yp) 


where  V„  and  Vp  are  the  band  parameters,  which  take  into 
account  the  variation  of  band  edge  with  doping  and  com¬ 
position,  the  shape  of  energy  band,  and  the  distribution 
statistics. 

Equations  (1)— (3)  can  be  generalized  and  expressed  in 
the  form 

V-{aVu)=f  (6) 

u  =  (i p,  v,  w)  (7) 

a  =  {e.  toe1**™,  fye1-***')  (8) 

/=  (n-p-N£+N;,^  +  R,ft  f  /?). 


Note  that  (6)-(9)  denote  three  individual  equations  in  the 
same  mathematical  form. 

If  the  Gummel  iteration  method  is  assumec'  mg 

the  Galerkin’s  analysis  of  (6),  the  finite  element  discreti¬ 
zation  results  in  the  formulation  [8] 


S  (K,J  +  G,j)Auj  =  -  L  KfjUj  +  f)  ( 10) 
i  -  ■  V  j  = 1  / 

where  the  stiffness  matrices  and  load  vector  are  defined  as 
*,=  55  aV^  •  V*,  ds  (11) 

(,2> 


/*  5  J/*.-*  (i3) 

and  <f>,  is  the  shape  function,  usually  chosen  as  a  low  order 
polynomial.  Note  that  there  is  no  derivative  term  of  a  in¬ 
volved  in  the  term,  however,  which  may  appear  in  cer¬ 
tain  cases  as  stated  in  Appendix  C.  Here,  the  derivative 
of  a  is  implicitly  included  in  the  decoupled  iteration 
scheme,  where  the  Poisson  equation  is  solved  first  by  as¬ 
suming  known  quasi-Fermi  potentials,  then  continuity 


equations  are  solved  with  the  potential  given  from  the  pro¬ 
ceeding  step.  And  the  Gummel  iteration  is  repeated  until 
self-consistent  values  of  all  unknowns  are  obtained. 

For  finite  element  analysis,  the  main  concern  is  the  ac¬ 
curacy  of  the  quadrature  in  the  assembly  of  stiffness  ma¬ 
trices  and  right-hand  side  vector  of  (10),  i.e.,  assembling 
(1 1 )— ( 13),  in  which  the  complex  calculations  are  involved 
in  determining  the  values  of  the  function  to  be  integrated. 

Conventionally  a  set  of  low  order  polynomials  [2]  are 
also  used  to  approximate  the  carrier  densities,  n  and  p, 
(e(^  +  |/")),  and  +  which  can  be  written  as 


n  =  2  n,  4>j 

ite 

(14) 

p  =  Z /?,<*>, 

ite 

(15) 

e(*  +  V„)  =  £  (gt* +  *'»))  <£. 

ite  1 

(16) 

07) 

ite 


where  index  i  is  the  alternation  of  the  vertices  of  the  tri¬ 
angular  element  e.  Let’s  take  an  example  of  assembling 
the  matrix  element  Gf/e  of  Poisson  equation  due  to  ele¬ 
ment  e. 

G*e = 55 (n +p)Mjds 

e 

where  two  methods  can  be  applied  for  the  approximation 
of  such  integral  functions.  On  the  one  hand,  the  direct 
integration  is  an  intuitional  choice;  on  the  other  hand,  the 
Gaussian  quadrature  rules  with  various  orders  are  appli¬ 
cable.  If  the  direct  integration  method  is  applied,  the  re¬ 
sultant  approximation  will  be  simply  in  the  form 

Gf/*  =  5  5  ^  (n  +  P)k4k4>i4>j  ds 

e 

-  S  (n  +  p)k  5  5  ds.  ( 19) 

e 

It  is  clearly  seen  from  the  above  that  using  low  order 
polynomials  to  approximate  the  carrier  densities,  n  and  p , 
(e(*  +  v’))t  ancj  -*  +  v?)}  wjn  no  doubt  result  in  large 
integration  errors  because  these  variables  inherently  ex¬ 
hibit  a  highly  nonlinear  variation  in  a  small  region,  lead¬ 
ing  to  an  inaccurate  solution  as  well  as  a  poor  conver¬ 
gence  rate,  except  that  the  mesh  is  intensively  refined. 

The  Gaussian  quadrature  requiring  a  least  number  of 
evaluations  to  achieve  a  maximum  accuracy  is  ideally 
suited  for  our  use  [8].  If  the  Gaussian  quadrature  rules  are 
applied,  we  might  need  values  of  n,  p ,  {e{i  +  v')).  and 
(e[  +  Vf))  at  certain  points  in  addition  to  those  at  vertices 

of  the  element  in  order  to  have  an  accurate  high  order 
approximation,  which  could  be  interpolated  from  the 
known  values  at  vertices.  However,  it  is  obvious  that 
using  the  same  low  order  interpolation  polynomials  like 
( 14)— ( 1 7)  in  Gaussian  quadrature  will  suffer  from  the  same 
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low  accuracy  problem  as  that  which  occurs  in  the  direct 
integration  method.  Thus  an  approximation  of  the  inte¬ 
grand  in  the  Gaussian  quadrature  based  on  a  nonlinear 
functional  fitting  would  be  more  desirable. 

III.  New  Finite-Element  Discretization  Employing 
the  S-G  Scheme 

In  order  to  take  advantage  of  the  simplicity  and  accu¬ 
racy  of  the  S-G  scheme  associated  with  the  finite  differ¬ 
ence  discretization  and  to  avoid  the  shortcomings  of  the 
classical  finite  element  integration  method,  we  employ  a 
seven-point  Gaussian  quadrature  rule  in  conjunction  with 
the  S-G  scheme,  while  assembling  the  finite-element  stiff¬ 
ness  matrices  and  right-hand  side  vector.  This  new  scheme 
can  be  outlined  as  follows. 

For  an  arbitrary  triangular  element,  e,  as  shown  in  Fig. 
1 ,  the  seven-point  Gaussian  quadrature  rule  with  error  es¬ 
timate  of  0(h 4)  [8]  reads 

/=  (  fds  =  2A  Z  (20) 

where  A  is  the  area  of  element  e;  Lx  and  Ln  are  natural 
coordinates  of  the  triangular  element;  the  seven  sampling 
points  on  an  element  are  taken  at  the  veitices  of  triangle, 
the  middle  points  of  the  edges  and  the  center  of  triangle, 
respectively,  as  illustrated  in  Fig.  1;  and  the  correspond¬ 
ing  weighting  coefficient  for  Gaussian  quadrature  at  the 
fcth  point,  wk,  is  listed  in  Table  I. 

If  the  quadrature  rule  (20)  is  applied  to  the  finite  ele¬ 
ment  integration  of  (11)-(13),  a  general  form  can  be  de¬ 
scribed  as 

j  j  G{n,p,  e^  +  ^\  e(-*  +  y’))ds 
< 

=  2A  Z  G(nk,  pk ,  (eu'v',))k,  (el~+  +  v',))k)  wk 

(21) 

where  nk,  pk,  (e{~*  +  Vp,)k,  ( k  =  /,  m,  n,  p) 

can  be  evaluated  using  a  S-G  formula-like  function.  For 
the  approximations  of  nk  and  pk,  we  follow  the  S-G  scheme 
and  assume  the  electric  field  and  the  current  densitv  be¬ 
tween  two  neighboring  grid  point  /  and  j  are  constants, 
which  yields  the  interpolation  formula  based  on  current 
density  equations  (see  Appendix  A  for  details). 

«(*)|*  =  l1  ~  «(**.  W<j)]n,  +  ^hj)n}, 

Xk([xi,Xj]  (22) 

p(*)\k  =  n  -  #(**■  -ty;j)]pi  +  g(xk,  -a  i^pj, 

xkt[x„  Xj)  (23) 

1  -ft! 

s(.r,  y)  =  - (24) 


i 


Fig.  1.  A  triangular  element  and  its  sampling  points  of  Gaussian  quadra¬ 
ture. 


TABLE  I 

The  Weightin'?  Coefficients  for  the  Seven-Point  Gaussian 
Quadrature 


>  i 

k 

1  tn 

0 

P 

2w, 

3/60  3/60 

3/60 

S/60  S/60 

1/60 

27/60 

a h  =  ij  +  Vnj  -  (*,  +  Vni)  (25) 

hy  Xj  Xr  ( 26 ) 

For  the  approximations  of  (^(^  +  l/"))  and  ( +  ), 

Bank  [9]  pointed  out  that  the  following  three  schemes  can 
be  applied: 

(1)  linear  interpolation  from  the  nodal  values 
(e(*  +  1'")),  and  (e(*  +  l'")); 

(2)  linear  interpolation  from  the  nodal  values  (\p  +  Vn)i 
and  +  Vn)j 

(3)  interpolation  by  using  Bernoulli  function 

and  suggested  that  scheme  (3)  is  more  accurate  than  the 
others  in  the  sense  of  nonlinear  functional  fitting  Based 
on  device  physics  considerations,  we  derive  an  interpo¬ 
lation  rule  similar  to  the  scheme  (3),  which  can  be  ex¬ 
pressed  as 

(e+  +  y')k  =  B{  A^)e<^K*>  (27) 

(e~*  +  y')k  =  +  (28) 

with  Bernoulli  function  defined  as 

*(*)  =  pr—j--  (29) 

Concerning  the  determination  of  the  function  value  at 
the  central  point,  we  take  an  average  of  the  three  inter¬ 
polated  values  which  are  obtained  by  applying  the  inter¬ 
polation  rule  described  above  along  the  respective  line 
connecting  a  vertex  and  the  middle  point  of  the  corre¬ 
sponding  edge  opposite  to  it,  i.e.,  i  —  n,j  —  /,  and  k  - 
m  in  Fig.  1 .  In  an  appropriately  refined  mesh  which  when 
used  in  the  new  scheme  is  always  more  coarse  than  that 
for  the  classic  finite  element  method,  the  differences 
among  these  three  interpolated  values  are  negligible  since 
we  have  a  more  accurate  functional  fitting  expression  cho¬ 
sen  for  interpolation. 
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For  clarity,  let’s  formulate  the  quadratures  for  Poisson 
equation  and  electron  continuity  equation  in  the  foHowm0 
(decoupled  iteration  method  assumed): 

(1)  Poisson  equation:  The  calculation  of  Kfe  is  sim¬ 
ple,  since  a  =  e  is  constant.  Here  only  formulations  for 
Gfj'e  and  ft e  are  given,  which  read 

Gf/e  =  j  j  («  +  p)<t>i<t>j  ds 

f 

7 

=  2A  ^  K  +  Pk)  ( <t>i)k  ( 4>J  )kwk  (30) 

ft'  =  ff  ~  P  ~  Ns)d>i  ds 

t 

7 

=  2A  E  K  -  Pk  -  Nsk)  (<t>,  )  wk  (31) 

k  *  I 

where  nk,  pk  (k  =  l,  m,  n,  p)  can  be  evaluated  from  (22) 
and  (23),  respectively,  and  (<£,  )*  =  <t>,  (xk,  yk),  ( 4>j)k  = 
<t>)  ( xk ,  yk). 

(2)  Continuity  equation: 

KY  =  j  j  P'eU  +  ^VtMjds 

e 

7 

=  2A  £  *i*(**  +  K,)t  (VA  )fc  (Vd>j \wk  (32) 


where  V<f>,  and  V<fy  are  constant  vectors,  if  linear  shape 
functions  are  chosen  for  triangular  element. 

For  simplicity,  let  us  assume  the  S-R-H  recombination 
dominates  the  recombination  mechanism 


f  =  R  = 


np  -  n, 


fn{p  +  «.•)  +  Tp{n  +  n,) 


(33) 


df  _  (P  +  nj)(jnP  + ;  _ (*  +  k„) 


7 


^  [r„(p  +  «,)  +  rp(n  +  /*,-)] 

=  R(n,p)nire(*  +  t'l') 
then  the  matrix  elements  are  given  by 

GY  R(n,  p)nir{e^  +  v’,))<t>i<t>J  ds 


(34) 


7 

=  2A  R(nk,  pk)ntr(e('i'*K))k  (<f>,)t  (<fy)tw* 

(35) 

fY  =  \\  R<t>,ds  =  2A  tZ  R(nk.  pk)  (*,.),*■*  ('<M 


where  for  fc  =  (/,  m,  n,  p),  n*.  pk.  and  (eu  *  K,>)t  are 
calculated  by  (22)-(29). 

In  the  sense  of  finite  element  analysis,  the  FEM-QSG 
formulation  conforms  to  the  Galerkin  method,  where  the 


shape  functions  serve  as  the  weighting  functions.  In  sit¬ 
uations  where  the  dominant  drift  term  is  strong  enough  to 
cause  spurious  oscillations,  we  can  easily  modify  the 
weighting  functions  in  the  FEM-QSG  formulation  to 
achieve  upwind  effect  without  any  inconsistency  [4],  (See 
Appendix  C  for  a  detailed  discussion.) 

In  order  to  take  advantage  of  the  simplicity  in  deriva¬ 
tion  and  its  consistency  with  the  finite  element  theory,  the 
definition  of  current  proposed  by  Barnes  et  al.  [2]  is 
adopted  in  our  current  calculation,  by  which  the  local  as 
well  as  global  current  conservations  are  preserved.  It  must 
be  stressed  that  in  the  FEM-QSG  formulation  the  modi¬ 
fication  on  the  interpolation  of  n,p,  e *  +  v\  for  the 

quadratures  would  not  perturb  the  current  conservation 
property  and  the  consistency  with  the  finite  element  the¬ 
ory. 

This  new  discretization  scheme  not  only  uses  an  accu¬ 
rate  seven-point  Gaussian  quadrature  rule  instead  of  the 
classical  integration  method  but  also  applies  the  S-G 
scheme  to  calculating  the  values  of  quadrature  points 
which  are  not  on  the  vertices  of  an  element  so  that  the 
approximation  of  carrier  densities,  (e(*  +  l'',))  and 
(e(~'i  +  l^))  yields  an  accurate  exponential  functional  fit¬ 
ting  form.  As  a  result,  the  convergence  rate  would  be  sig¬ 
nificantly  enhanced  and  the  accuracy  of  solution  would  be 
obviously  improved  due  to  the  low  finite  element  discre¬ 
tization  error  generated  by  the  new  method. 

IV.  Simulation  Results  and  Discussions 

In  the  BIM  discretization,  an  equivalent  central  finite 
difference  discretization  scheme  is  usually  applied  to  the 
equation  assembly  of  Poisson  equation  in  a  two-dimen¬ 
sional  domain,  which  results  in  a  solution  with  compara¬ 
ble  accuracy  to  that  generated  by  the  classical  finite  ele¬ 
ment  method  using  piecewise  linear  finite  elements  [10]. 
Additionally,  the  linear  interpolation  on  the  source  terms 
(the  right-hand  side  terms)  of  Poisson  equation  by  the  BIM 
will  inevitably  introduce  a  larger  discretization  error  than 
that  by  the  S-G  interpolation  in  our  method,  especially  for 
the  strong  coupling  cases.  On  the  solution  of  current  con¬ 
tinuity  equations,  a  hybrid  finite-difference  finite-element 
technique  in  conjunction  with  the  S-G  scheme  is  applied 
in  the  BIM  scheme.  From  the  viewpoint  of  the  general¬ 
ized  finite-element  method  [1 1],  it  is  found  that  there  is  a 
natural  correspondence  between  the  BIM  scheme  and  our 
new  method  on  the  discretization  of  continuity  equations. 
A  more  detailed  comparison  is  presented  in  Appendix  B. 

To  demonstrate  the  effectiveness  and  accuracy  of  this 
new  discretization  method  (method  II)  compared  to  the 
classical  one  (method  I),  we  give  three  examples.  From 
the  theoretical  analysis,  we  gather  that  the  convergence 
rate  is  strongly  dependent  on  the  accuracy  of  discretiza¬ 
tion.  We  thus  consider  that  a  detailed  comparison  with  the 
classical  finite  element  can  show  us  to  what  extent  that 
the  discretization  errors  have  influence  on  the  conver¬ 
gence  rate. 

First  of  all,  a  simple  abrupt  p-n  junction  is  simulated. 
Two  meshes  used  for  simulations  are  illustrated  in  Fig.  2. 
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Fig.  2.  Meshes  used  for  simulations  of  an  abrupt  p-n  junction  with  grid- 
points  of  (a)  NP  =  37,  (b)  NP  =  640. 


The  coarse  mesh  is  used  as  a  reference  for  the  following 
evaluations.  The  mesh  with  fine  grid  is  a  typical  mesh 
example  for  obtaining  the  optimal  solution.1 

Fig.  3  shows  the  number  of  iterations  required  by 
method  II  to  reach  convergence  and  that  required  by 
method  I,  as  a  function  of  impurity  concentration.  Fig.  4 
compares  the  minimum  number  of  mesh  points  required 
for  convergence  by  method  II  with  that  by  method  I,  also 
as  a  function  of  impurity  concentration.  As  seen  in  Figs. 
3  and  4  the  results  indicate  that  the  integrations  by  method 
II  are  so  accurate  that  the  convergence  speed  is  almost 
independent  of  the  impurity  density  and  the  grid  size  of 
the  mesh.  On  the  contrary,  the  convergence  rate  of  the 
simulation  discretized  by  method  I  is  very  sensitive  to  the 
impurity  density  and  the  fineness  and  the  layout  of  the 
mesh  due  to  large  integration  errors  involved  in  the  ap¬ 
proximation  of  nonlinear  functions  during  assembling  the 
element  matrices,  especially  for  highly  doped  regions. 

Fig.  5  shows  the  decreasing  rate  of  the  error  function 
of  Poisson  equation  (at  zero  bias).  It  is  obvious  that  the 
convergence  rate  is  quadratic  for  method  II.  Besides,  it 
does  not  need  dt  mping  scheme  in  Newton  iteration  except 
the  first  few  iterations,  whereas  the  convergence  rate  is 
almost  a  linear  one  for  method  I. 

In  Fig.  6,  we  compare  the  discretization  error,  ERR,2 
of  Poisson  equation  by  method  II  with  that  by  method  I, 


'Here  the  optimal  solution  is  obtained  by  keeping  reducing  the  grid  size 
by  half  each  time  (i.e.,  refining  the  mesh)  until  the  relative  deviation  be¬ 
tween  solutions  of  two  successive  discretizations  is  less  than  1  percent. 

JThe  discretization  error  ERR  is  defined  as  the  relative  error  between 
the  exact  solution  and  numerical  solution. 


Fig.  3.  The  number  of  iteration  NUM  required  for  convergence  versus  im¬ 
purity  concentration  A/D  (N,,  =  1  x  iO15  cm*3,  zero  bias). 


Fig.  4.  The  minimum  number  of  mesh  points  Np„,„  required  for  conver¬ 
gence  versus  impurity  concentration  N„  (NA  =  1  x  10”  cm'3,  zero 
bias). 


Fig.  5.  The  decrease  of  the  error  function  of  Poisson's  equation  versus  the 
number  of  iterations  (NP  =  37,  zero  bias)  by  (1)  method  I  and  (2)  method 
II. 


as  a  function  of  minimum  mesh  size.  And  we  show  in  Fig. 
7  the  discretization  errors  of  continuity  equations  of  elec¬ 
trons  and  holes  by  method  II  with  that  by  method  I.  again 
as  a  function  of  minimum  mesh  size. 

Figs.  6  and  7  confirm  that  by  method  II  we  can  obtain 
a  more  accurate  solution  than  that  by  method  I  on  the  same 
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Fig.  6.  The  discretization  errors  ERR  between  exact  solution  and  numer¬ 
ical  solution  of  Poisson’s  equation  as  a  function  of  minimum  mesh  size 
A Amin  ( Na  =  1  x  10”  cm  ,  Nd  =  1  x  10”  cm"3,  zero  bias). 
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Fig.  7.  The  discretization  errors  of  continuity  equations  of  electron  and 
hole  between  exact  solution  and  numerical  solution  as  a  function  of  Axm,„ 
(/V„  =  1  x  10”  cm"3,  Na  =  5  x  10”  cm"3,  Klpp  =  0.6  V). 

grid,  in  particular  for  the  Poisson  equation  case.  In  addi¬ 
tion,  the  results  illustrate  that  applying  method  II  to  con¬ 
tinuity  equations  will  allow  us  to  take  the  same  advan¬ 
tages  of  the  S-G  scheme  as  those  well  known  in  the  finite 
difference  counterpart.  Therefore,  this  new  finite  element 
scheme  is  competitive  with  the  finite  difference  method 
with  exponential  functional  fitting  in  the  sense  of  discre¬ 
tization  accuracy. 

Fig.  8  shows  the  simulation  potential  distribution  and 
electron  and  hole  concentration  plots  of  the  abrupt  p-n 
junction  of  the  optimal  solution1  compared  with  those  ob¬ 
tained  by  methods  I  and  II  using  the  mesh  shown  in  Fig. 
2(a).  We  can  see  that  the  solution  calculated  by  method 

II  is  closer  to  the  optimal  solution  than  that  by  method  I 
on  the  basis  of  same  simulation  mesh. 

In  order  to  test  the  effectiveness  and  efficiency  of  the 
new  method  for  complex  device  cases,  we  further  simu¬ 
lated  a  typical  bipolar  transistor  and  a  MOS  transistor. 

The  simulated  bipolar  transistor,  whose  geometric 
structure  is  shown  in  Fig.  9(a),  is  doped  with  an  impurity 
concentration  profile  as  plotted  in  Fig.  9(b).  Tables  II  and 

III  list  the  relative  deviation,  RD3  of  the  calculated  results 

3The  relative  deviation,  RD.  is  defined  as  the  maximum  relative  devia¬ 
tion  of  the  calculated  results  from  the  optimal  solution. 


Fig.  8.  Plots  of  (a)  the  potential  distribution,  (b)  the  electron  concentra¬ 
tion,  and  (c)  the  hole  concentration  of  along  the  central  line  of  the  abrupt 
p-n  junction  of  (1)  the  optimal  solution,  solutions  obtained.  (2)  by  method 
1,  and  (3)  by  method  II. 

v  ith  respect  to  the  optimal  solution  for  the  bipolar  tran¬ 
sistor  biased  at  Vc  =  1.6  V  and  Vg  =  0.9  V  with  various 
meshes  by  methods  I  and  II.  The  results  manifest  that  by 
using  method  II  the  solution  has  higher  accuracy  than  that 
by  method  I.  In  other  words,  under  the  same  accuracy 
requirement,  method  II  can  use  a  coarser  mesh  than  that 
used  by  method  I  to  reach  the  same  goal. 
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(a) 


Fig.  9.  (a)  The  device  structure  and  (b)  the  doping  profile  along  the  line 
of  y  =  0  of  the  simulated  bipolar  transistor. 


TABLE  II 

The  Relative  Deviation  of  the  Calculated  Results  from  the 
Optimal  Solution  for  the  Bipolar  Transistor  (at  Vc  =  1.6  V.  Vb  = 
0.9  V)  using  Various  Meshes  (Method  I) 


NP  |  554 

90* 

RD(%)v 

2.775 

2.774 

RD(%)  n 

8.722 

8.720 

RDOMp 

3-512 

5.731 

TABLE  III 

The  Relative  Deviation  of  the  Calculated  Results  from  the 
Optimal  Solution  for  the  Bipolar  Transistor  (at  K,  =  1.6  V.  V„  = 
0.9  V)  using  Various  Meshes  (Method  II) 


Secondly,  a  MOS  transistor,  whose  structure  is  shown 
in  Fig.  10,  is  simulated  as  well.  Fig.  11  shows  the  con¬ 
vergence  rate  of  Poisson  equation  for  MOS.  We  found  out 
that  method  II  also  held  a  quadratic  convergence  rate.  We 


(„m)x 


(a) 


Fig.  10.  (a)  The  device  structure  of  the  simulated  MOS  transistor  and  (b) 
the  mesh  used  for  simulation. 


Fig.  1 1.  The  decrease  of  the  error  function  of  Poisson's  equation  versus 
the  number  of  iterations  of  Poisson's  equation  for  the  MOS  (zero  bias) 
by  (I)  method  I  and  (2)  method  II. 


tabulate  the  RD  of  the  MOS  example  by  methods  I  and  II 
in  Tables  IV  and  V,  respectively,  and  the  results  also  re¬ 
veal  that  the  accuracy  of  sotution  by  method  II  is  higher 
than  that  by  method  I.  Fig.  12  shows  perspective  plots  of 
electron  and  hole  distribution  for  the  BJT  (at  Vc  =  1.6  V, 
Fft  =  0.9  V). 

We  can  easily  find  out  from  the  examples  demonstrated 
above  that  the  advantages  of  the  new  discretization  method 
applied  equally  well  to  a  simple  abrupt  p-n  junction  and 
the  complicated  device  structures. 
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Fig.  12.  The  perspective  plots  of  electron  and  hole  distributions  of  the 
simulated  BJT  (  Vc  =  1 .6  V,  V„  =  0.9  V ). 

TABLE  IV 

The  Relative  Deviations  of  the  Solutions  for  n-MOS  (at  V,  =  1.0 
V,  Vd  =  0.5  V )  Estimated  Along  the  Lines  of  x  =  0.04  and  0.06 
(Method  II) 


NP 

220 

375 

RD(%) 

*•0.04 

1.28 

1.08 

V 

*•0.06 

146 

1.60 

RD  (%) 

*^><X 

3.84 

4.73 

n 

*-0  06 

3  86 

4.70 

TABLE  V 

The  Relative  Deviations  of  the  Solutions  for  n-MOS  (at  f',  =  1.0 
V,  Vj  -  0.5  V)  Estimated  Along  the  Lines  of  x  =  0.04  and  0.06 
(Method  I) 


NP 

|  340 

295 

492 

RD<%) 

*■0.04  403 

3.04 

3.71 

V 

*■0  06  4  68 

422 

402 

CL 

RD(%) 

I-O.04  1.6  S 

8  33 

596 

n 

*-0  06  7.80 

766 

5  19 

_ 
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The  weakness  of  the  new  method  stems  from  the  con¬ 
stant  derivative  assumption  of  i p,  4>„,  ■  ■  ■  imposed  on  the 
S-G  scheme.  The  approximation  is  poor  by  the  S-G 
scheme  at  the  transition  layer  where  a  large  variation  of 
electrical  field  or  current  density  occurs.  The  remedy  to 
this  is  refining  the  mesh  properly  by  self-adaptive  algo¬ 
rithm. 

V.  Conclusions 

A  simple  and  accurate  finite  element  discretization 
method  has  been  developed  in  the  two-dimensional  semi¬ 
conductor  device  analysis  program  SDA-1.  The  key  of 
this  method  lies  in  (1)  accurate  interpolation  rules,  which 
are  derived  based  on  simple  device  physics  considera¬ 
tions,  are  applied  to  the  functional  approximations  of  n, 
p  +  .  .  .  >  necessary  in  the  finite  element  assembly 

of  stiffness  matrices  as  well  as  source  terms  of  Poisson 
equation  and  continuity  equations  in  a  consistent  way 
rather  than  just  concentrating  on  the  discretization  of  flow 
term,  V  •  (aV«),  of  continuity  equation  like  the  BIM-SG 
scheme  does.  (2)  The  inherent  simplicity  and  flexibility 
in  the  finite  element  formulation  make  the  new  method 
applicable  to  multidimensional  problems.  (3)  The  sim¬ 
plicity  of  embedding  the  S-G  scheme  in  the  quadrature  of 
finite  element  assembly  lends  itself  to  all  kinds  of  finite 
element  methods  employing  various  elements,  shape 
functions  and  weightings.  The  resultant  exponential  func¬ 
tional  fitting  by  the  new  method  avoids  high  discretization 
errors  usually  incurred  by  the  classical  finite  element  dis¬ 
cretization  method.  Using  this  new  method,  we  not  only 
obtain  solution  with  high  accuracy  but  also  speed  up  the 
convergence  rate  significantly.  Moreover,  as  the  theoret¬ 
ical  prediction,  the  numerical  experiments  confirm  that  the 
convergence  rate  strongly  relies  on  a  proper  function'd  fit¬ 
ting  and  an  accurate  numerical  integration  in  the  appli¬ 
cation  of  the  finite  element  method  to  the  semiconductor 
equations. 

Appendix  A 

The  Application  of  the  S-G  Scheme  to  Finite 
Element  Method  Based  on  Device  Physics 
Considerations 

Considering  the  current  density  equation  for  electrons, 
in  terms  of  different  dependent  variables,  they  can  be  ex¬ 
pressed.  respectively,  as 


Jn  =  ~  Vn) 

(Al) 

Jn  = 

(A2) 

According  to  the  Schafetter-Gummel  (S-G)  scheme,  if 
the  derivatives  of  \p,  v,  •  ■  •  ,  between  point  /  and  j  are 
assumed  constant,  the  S-G  approximation  of  the  current 
density  equation  projected  on  i  -  j  line  by  central  differ¬ 
ences  reads 

Jn  =  h  —  (A3) 


476 


IEEE  TRANSACTIONS  ON  COMPUTER-AIDED  DESIGN,  VOL.  S.  NO.  5.  MAY  1989 


where  B{ x)  =  x/(e*  —  1 )  is  the  Bernoulli  function;  Ji  is 
the  average  value  of  between  point  i  and  j  [9].  It  seems 
difficult  to  derive  a  similar  expression  like  (A3)  by  the 
finite  element  discretization.  However,  the  similarity  be¬ 
tween  (A3)  and  (A2)  suggests  that  e^B{  A\pj, )  is  a  proper 
interpolation  formula  for  evaluating  e*  at  certain  point  be¬ 
tween  points  i  and  j\  or  in  a  more  exact  term,  the  e *  value 
at  midpoint  M  between  points  i  and  j  can  be  approximated 
by 

(e\=  e^B(A^).  (A4) 

The  carrier  densities,  n  and  p,  are  well  known  for  hav¬ 
ing  an  exponential  functional  dependence  on  \p .  The  lin¬ 
ear  interpolation  will  obviously  result  in  an  inaccurate  es¬ 
timation.  However,  we  can  further  extend  the  idea  behind 
the  S-G  scheme  to  the  evaluation  of  carrier  concentration, 
at  point  k  located  in  between  two  neighboring  gridpoints 
i  and  j,  in  terms  of  the  given  carrier  concentrations  n,  and 
rij  at  these  two  terminal  point  as  follows. 

Based  on  the  constant  derivative  assumption,  the  cur¬ 
rent  density  equation  (Al)  in  a  one-dimensional  case  re¬ 
duces  to  a  two-point  boundary  value  problem: 


~  +  Etna  =  — 

dx  M/i 

(A5) 

with  the  boundary  conditions 

n(xi)  =  n, 

( A6) 

n(.Xj)  =  rij 

(A7) 

where  £cfr  =  - 

The  general  solution  n  can  be  expressed  as 

n  =  n,  ( 1  -  g)  +  Tijg 

1 

hv 

1  —  ey 

(A8) 

(A9) 

where  htJ  =  x,  -  x,,  y  =  [ ~  'Z',  ].  Then  the  equivalent 
S-G  formula,  (22),  can  be  obtained  and  (23)  can  be  de¬ 
rived  in  the  same  way. 


Appendix  B 

The  Comparison  of  the  New  Method  with  the 
BIM-SG  Scheme 

A  triangular  (or  rectangular)  domain  partition  process 
is  usually  employed  in  the  BIM  scheme  first.  Then  boxes 
are  constructed  around  every  gridpoint  by  taking  the  mid- 
perpendicules  of  line  segments  directed  to  its  neighboring 
points.  If  the  Green’s  theorem  is  applied  to  (6)  on  a  box, 
we  have 

|  oVu  dl  =  j  f  ds  (Bi ) 

or 

2  aVu  ■  n,l,  =  fS  (B2) 


where  /i,  is  the  unit  normal  vector  of  box  edge  l,  is  the 
length  of  the  edge;  S  is  area  of  the  box. 

If  the  BIM  scheme  is  applied  to  Poisson  equation  on  the 
example  mesh  shown  in  Fig.  13,  a  standard  central  dif¬ 
ference  expression  would  be  obtained: 

a(u,  +  u2  +  m3  +  k4  -  4 u0)  =  fuh-.  (B3) 

This  shows  an  equivalence  of  the  BIM  scheme  on  Poisson 
equation  to  the  standard  central  difference  with  accuracy 
of  second  order.  From  other  point  of  view,  the  application 
of  classical  FEM  to  Poisson  equation  on  the  same  mesh 
will  yield  a  same  central  difference  expression.  Thus  the 
BIM  scheme  is  equivalent  to  the  classical  FEM  on  the 
discretization  of  Poisson  equation.  Although  we  use  a 
rectangular  mesh  here,  the  derivation  holds  generally  true 
on  other  mesh  geometries  [10], 

In  Section  III,  we  have  shown  the  superiority  of  the  new 
method  to  the  classical  FEM  due  to  the  accurate  assembly 
of  the  stiffness  matrices.  Besides,  we  introduce  the  S-G 
scheme  to  the  Jacobian  and  right-hand  side  vector  via  J  j 
(n  +  p )  ds  and  J  j  (n  —  p  —  Ns)  ds.  Compared 
to  the  BIM  scheme  where  the  source  term  is  approximated 
simply  by  f0  times  the  box  area,  the  more  ac  jrate  treat¬ 
ment  of  the  right-hand  side  vector  in  the  new  method  will 
definitely  further  reduce  he  discretization  error.  Thus  the 
new  method  does  not  need  mesh  as  fine  as  that  usually 
used  by  the  BIM  sche  ,e  to  achieve  a  same  convergence 
rate  in  solving  Poisson  equation. 

As  for  continuity  equation,  the  BIM  scheme  incorpo¬ 
rates  the  S-G  current  formula  to  improve  the  discretiza¬ 
tion  accuracy.  For  convenience,  in  the  following  compar¬ 
ison  we  consider  the  contribution  to  the  Jacobian  only 
from  the  triangule  A015  as  indicated  in  Fig.  13  by  the 
BIM-SG  scheme,  which  gives 


J\  '  l\  +  J2  '  lz  —  M01  e't°^(A'/'oi ) 


V,  -  v0 


+  *>5  e*°BMw) 


"5  ~ 


=  (B4; 

where  for  the  mesh  in  Fig.  13,  /,  -  h/2  and  I2  -  0.  (F01 
other  cases,  /(,  l2  may  not  be  zero;  however,  the  compar 
ison  still  holds  valid.) 

If  the  new  method  is  applied,  according  to  (32).  th< 
discretization  of  continuity  equation  on  element  A01! 
reads 

2  K^vj  =  2AJi  2  wk(e*)  ■  2  (V4>,\  (V*,) 

/=0.  1.5  k  j  K 


Assume  the  shape  function  is  linear,  then 


V#,  •  V<*  = 


b,bj  +  c,Cj 
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Fig.  13.  An  example  rectangular  mesh. 


where  A  is  area  of  element  A015,  and  (B5)  becomes 


;  =  0, 1,5 


.  s  ,  KyV}  =  a  2  2W*(e*)4  s 


(bibj  +  c,c; 


V  4A 


(B7) 


where  is  average  value  of  p.  Moreover,  from  the  linear 
shape  function,  we  can  easily  prove  that 


2  b,b‘  +  C,Cj 
j  4A 


v0  -  v, 


(B8) 


Therefore,  (B7)  can  be  further  expressed  in  the  form 

where 


Z>2wk(e\ 


E 

Jfc  =  0. 1.5 


■  5(a^,5)  +  i  (^v  (B10) 


Equation  (BIO)  is  equivalent  to  a  weighting  scheme  over 
( e *)  with  weights  on 

(1)  three  vertex  values  ( e *)  by  3/60; 

(2)  midpoint  values  of  three  edges  (corresponding  to 
(e^)0  B{  A\f/0l )  in  the  S-G  scheme)  by  8/60; 

(3)  the  value  at  the  central  point  (e^)p,  which  is  also 
an  average  of  the  midpoint  values  and  the  vertex 
values,  by  27/60. 

Comparing  (B9)  with  (B4),  we  will  find  that  the  idea 
behind  the  S-G  scheme  is  embedded  in  the  new  method 
in  a  different  form  wher°  the  weighting  average  of  (e*) 
is  taken  on  the  element  rather  than  on  the  edges  by  the 
BIM-SG  scheme.  Thus,  the  new  method  is  at  least  as  good 
as  the  BIM-SG  scheme  in  the  sense  of  applying  the  S-G 
scheme  to  continuity  equation.  It  might  have  a  lower  dis¬ 
cretization  error  than  the  BIM-SG  scheme  does,  if  we 
count  the  factors  of  using  the  accurate  seven-point  Gauss¬ 


ian  quadrature  rule  and  weighting  the  average  of  a  on  ele¬ 
ment  instead  of  on  edges. 

Again,  the  new  method  has  a  more  accurate  treatment 
on  the  source  term  /  than  the  BIM-SG  scheme  does  by 
further  introducing  the  S-G  formula  to  the  terms,  j  j 
{df/du)  4>i <frj  ds  and  j  f/0,  ds,  where  the  BIM-SG  scheme 
simply  approximate  these  terms  by  linear  interpolations. 
In  this  respect,  if  decoupled  solution  method  is  used  and 
the  generation-recombination  terms  are  neglected  where 
the  continuity  equation  is  linear,  the  new  method  is  iden¬ 
tical  to  the  BIM-SG  scheme.  However,  if  the  generation- 
recombination  terms  are  taken  into  account,  then  the 
equation  being  a  nonlinear  one,  the  new  method  certainly 
would  have  a  better  performance  than  the  BIM  scheme 
due  to  a  lower  discretization  error  from  source  term. 

Appendix  C 

Relevant  Characteristics  to  the  Upwind  Effect 

If  a  in  (6)  is  a  function  of  u,  (10)  will  involve  a  first- 
order  derivative  term  (the  drift  term),  via 

V  •  [a(u)Vu]  =  a(u)V2u  +  Va(u)  ■  Vu  =  f  (Cl) 

The  drift  term  has  to  be  taken  into  consideration  at  the 
following  situations. 

(1)  A  coupled  solution  method  is  used. 

(2)  A  decoupled  iteration  method  is  used,  while  the 
mobility,  ju,  is  a  function  of  the  Fermi-potential  or  carrier 
density. 

As  is  well  known,  the  solutions  to  the  convective  u<ms- 
port  equation  like  (Cl)  by  the  standard  finite  difference 
method  or  the  classical  finite  element  method  are  often 
spoiled  by  spurious  oscillations.  To  preclude  such  oscil¬ 
lation  various  schemes  have  been  proposed  over  the  years 
for  finite  difference  discretization  method  and  finite  ele¬ 
ment  discretization  method.  In  the  finite  element  regime, 
artificial  diffusion,  quadrature,  and  the  Petrov-Galerkin 
method  are  common  techniques  utilized  to  achieve  the  up¬ 
wind  effect. 

The  application  of  the  S-G  scheme  has  been  shown  to 
be  equivalent  to  adding  certain  “upwind  term"  to  the 
standard  discretization  (5J.  In  the  new  method,  the  S-G 
scheme  is  incorporated  in  the  quadrature  and  an  accurate 
quadrature  rule  is  employed.  These  techniques  will  come 
up  to  enhance  the  upwind  effects  in  the  solution. 

Recently,  a  streamline  upwind/Petrov-Galerkin  method 

[4]  has  been  shown  effective  for  two-dimensional  convec¬ 
tion  dominated  flow  problems  by  choosing  weighting 
function  to  be  the  sum  of  the  standard  shape  function  and 
a  streamline  upwind  perturbation.  Since  the  new  finite 
element  formulation  FEM-QSG  conforms  to  the  standard 
Galerkin  method,  it  can  be  easily  converted  to  the  Petrov- 
Galerkin  formulation  by  simply  modifying  the  weighting 
function  vv,  based  on  the  shape  "unction  <f>,  without  causing 
any  inconsistency.  Therefore,  the  new  finite  element  dis¬ 
cretization  scheme  is  able  to  enhance  the  upwind  effects 
by  all  kinds  of  techniques. 
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ADIC-2.C:  A  General-Purpose  Optimization 
Program  Suitable  for  Integrated  Circuit 
Design  Applications  Using  the  Pseudo 
Objective  Function  Substitution 
Method  (POSM) 

GEN-LIN  TAN,  SHAO-WEI  PAN,  WALTER  H.  KU,  and  AN-JUI  SHEY 


Abstract— In  this  paper,  a  new  unconstrained  optimization  algo¬ 
rithm— POSM— is  presented.  This  algorithm,  requiring  neither  deriv¬ 
ative  calculation  nor  linear  search  step,  substitutes  the  objective  func¬ 
tion  by  a  second-order  approximated  formulation,  which  enhances  the 
convergence  rate  substantially  and  has  been  implemented  in  the  gen¬ 
eral-purpose  Analysis  and  Design  program  for  Integrated  Circuit, 
ADIC-2.C.  Its  high  efficiency  has  been  demonstrated  by  numerical  ex¬ 
amples  as  well  as  integrated  circuit  designs. 

I.  Introduction 

HE  COMMON  features  found  in  the  optimization 
problems  of  circuit  and  system  designs  are  a  limited 
number  of  parameters  to  be  optimized  and  a  huge  com¬ 
putation  load  that  is  incurred  due  to  an  excessively  com¬ 
plicated  objective  function  evaluation  need  (especially  for 
the  optimization  of  transient  characteristics).  For  the  ex¬ 
isting  optimization  algorithm,  calculation  of  the  deriva¬ 
tives  of  the  objective  function  or  its  equivalences  is  usu¬ 
ally  a  necessity.  Because  of  an  increasing  complexity  of 
the  circuit  and  the  high  nonlinearity  of  active  elements 
involved,  this  is  a  complicated  process  concerning  mod¬ 
ern  integrated  circuit  designs.  Thus  the  reduction  in  the 
number  of  function  evaluations  is  crucial  to  its  efficiency. 
Furthermore,  the  linear  search  is  an  essential  step  to  en¬ 
sure  the  convergence  of  Newton’s  iteration.  The  tech¬ 
niques  employed  during  the  linear  search  step  in  existing 
optimization  schemes  are  very  time  consuming.  There¬ 
fore,  it  is  our  intention  to  develop  a  new  algorithm  with 
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the  following  characteristics: 

(1)  no  need  for  derivative  calculation, 

(2)  reducing  the  need  of  linear  search  to  a  minimum, 

(3)  converging  in  fewest  iterations. 

Chen  [1]  has  developed  an  effective  optimization  al¬ 
gorithm  suitable  for  integrated  circuit  design  with  the  fol¬ 
lowing  characteristics:  (1)  making  full  use  of  the  evalu¬ 
ated  objective  function  values  (OFV)  in  much  the  same 
way  that  the  simplex  methods  [2]  do  and  (2)  estimating  a 
falling  direction  from  current  n  +  1  OFV’s  and  finding  a 
new  point  with  an  improved  objective  function  along  the 
falling  direction  in  a  manner  similar  to  the  Newton's 
method,  rather  than  blindly  taking  the  reflection  point  with 
respect  to  the  center  in  the  simplex  methods  (SM).  Thanks 
to  these  features,  Chen’s  method  proved  to  be  more  ef¬ 
fective  than  other  known  optimization  methods  [1].  How¬ 
ever,  the  resort  to  linear  search  is  still  frequent  in  Chen’s 
method  and  the  amount  of  iterations  is  huge.  Thus  we 
proposed  a  new  algorithm—POSM.  The  principal  idea  lies 
in  finding  the  minimum  point  of  the  pseudo  objective 
functions  which  are  formulated  to  approximate  the  exact 
objective  function  in  a  much  simpler  form  without  in¬ 
volving  complicated  calculation  instead  of  looking  for  the 
minimum  point  of  the  exact  objective  function  by  a  linear 
search  along  a  falling  direction.  Note  that  the  error  func¬ 
tion  can  be  approximated  by  a  first-,  second-,  or  even 
higher  order  polynomial  of  which  the  minimum  point  is 
easy  to  find,  and  a  series  of  optimization  methods  with 
various  degrees  of  complexity  can  be  developed  accord¬ 
ingly  following  this  idea.  Moreover,  the  modified  simplex 
method  in  conjunction  with  a  quadratic  interpolation 
scheme  is  incorporated  to  further  alleviate  the  necessity 
of  linear  search  step.  The  total  combination  leads  to  the 
pseudo  objective  function  substitution  method  (POSM). 

In  Section  II,  the  basic  concept  and  a  detailed  algo¬ 
rithmic  procedure  of  the  pseudo  objective  function  sub¬ 
stitution  method  are  described.  The  evaluation  results  of 
the  proposed  algorithm  by  a  few  typical  optimization 
problems  selected  from  literature  are  discussed  and  com- 
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pared  with  oth“r  well-known  algorithms  in  Section  III. 
Section  IV  presents  some  circuit  design  examples  illus¬ 
trating  the  great  suitability  of  the  application  of  POSM  in 
ADIC  for  integrated  circuit  and  system  design.  The  final 
section  presents  summary  and  conclusions  about  this  al¬ 
gorithm  and  comments  on  its  further  applicability. 


II.  Algorithm 


Generally  speaking,  the  circuit  optimization  can  be  per¬ 
formed  by  successively  modifying  the  parameter  vector  x 
=  (X[,  x2,  •  •  •  ,  x„)  until  the  calculated  response  u(x,  t) 
comes  close  to  the  desired  response  r(r)  within  a  certain 
tolerance  in  the  sense  of  a  least  pth  approximation  [3],  in 
short,  finding  an  x  such  that  the  objective  function  S(x) 
being  a  minimum 

S{x)  =  £  [ej(x)]P  (1) 

;*  i 


where  p  is  an  even  positive  integer,  m  is  the  number  of 
sampling  points,  and  e;(x)  is  the  corresponding  error 
function  at  the  y'th  sampling  point,  and  is  defined  as 

ej{x)  =  u(x,  t})  -  r{tj),  j  =  1,  2,  •  •  •  ,  m 


where  t,  for  example,  can  be  input  voltage  in  dc  charac¬ 
terization,  frequency  in  ac  analysis,  or  time  in  TRANsient 
evaluation;  u(x,  tj)  and  r(tj)  are  the  calculated  response 
with  parameter  vector  x  and  the  desired  response  at  the 
yth  sampling  point,  respectively.  The  basic  assumption  of 
POSM  is  that  es  can  be  approximated1  by 


ej  =  e“  +  S  btJ(x,  -  or")  +  £  Cy(x,  -  x")\ 


i  -  1 


j  =  1,  2,  •  •  •  ,  m 


(2) 


where  the  notation  e;  is  for  ej(x),  ef  for  ej(xM),  and  by 
and  Cy  are  approximation  coefficients  to  be  determined.  In 
order  to  determine  these  approximation  coefficients  a 
super-polygon  with  2n  +  1  vertices  in  n-dimensional  pa¬ 
rameter  space  is  constructed,  and  the  error  function  ek  = 
e;(x*),  and  the  objective  function  Sk  are  evaluated  at  these 
vertices  (Jfc  =  1,  2,  •  •  •  ,  2n  +  1 ,/  =  1,  2,  *  •  •  ,  m), 
from  which  xM  and  xmax,  the  parameter  vectors  of  the  best 
and  worst  point,  i.e.,  SM  (Smax)  with  the  smallest  (big¬ 
gest)  value,  among  the  2 n  +  1  points  can  be  selected  as 
a  reference  for  updating  the  super-polygon.  Putting  these 
values  (ek  and  ey)  into  (2),  we  can  formulate  a  set  of  2 n 
equations  in  terms  of  coefficients  b,y  and  c,y  (for  each  j ): 


*  =  1,  2,  •  •  •  ,  2 n  (3) 


'Here  the  mixed  terms  in  the  Taylor  series  expansion  are  neglected  for 
reducing  die  computation  load  when  calculating  those  coetticients.  Never¬ 
theless,  the  experimental  results  do  not  snow  high  sensitivity  to  those  mixed 
terms. 


where  wk  =  (Sk)  x/p  is  the  weight  of  the  fcth  equation 
[1],  Equation  (3)  can  be  expressed  in  a  matrix  form 


"(*l  r)  - 


(4) 


where  £  =  [wk(ek  -  ef  )]2nxm,  X  =  [x‘  -  x«]2nXn, 
Y=  [(x*  -  xf,)2]2nXn,B  =  (by]nxm,  and  C  =  [Cy]nXm. 
Ej,  Bj,  and  C;  denote  column  vectors  of  £,  B,  and  C,  re¬ 
spectively.  From  (4),  Bj  and  Cj  car.  be  solved  by 


=  [w(X|T)]+£;,  y  =  1,  2,  •  •  •  ,  m  (5) 

where  the  notation  [  ]+  represents  a  generalized-inverse 
operation  [4].  Note  that  we  only  need  to  calculate  B  and 
C  from  (5)  with  generalized-inverse  operation  once,  since 
w(Xj  Y)  in  (5)  is  independent  of  j. 

Based  on  the  second-order  error  function  approxima¬ 
tion.  a  pseudo  objective  function  $  corresponding  to  (1) 
can  be  formulated  in  terms  of  the  Taylor  series  expansion 
of  S(x )  at  SM  [5]  (see  Appendix  4). 

5(x"  +  Ax)  =  SM  +  PAxTBR 

+  0.5 p(p  -  1  )AxtBDBtAx  +  pYTCR 

+  0.5 p{p  -  1  )YtCDCtY 

+  p{p  -  l)AxrBDCrT  (6) 


where 

(e?r\ 

D  —  diag  [(etf)P_2,  •  •  •  ,  (^)P"2] 

Y  =  [Ax?,  Ax2,  •  •  •  ,  Ax?],  B  =  [bij]nxm, 


The  minimum  point  of  S  can  be  found  by  equating  the 
derivative  of  §  with  respect  to  Ax  to  zero 

«->  ■  -  »■  <7> 


From  (6)  we  have 

g(  Ax)  =  pBR  +  p(p  —  \)BDBtAx  +  2  pZCR 

+  p(p  -  1  )BDCtY 

+  2 p(p  -  1  )ZCDBT Ax 

+  2p(p  -  \)ZCDCtY  =  0  (8) 

where  Z  =  diag  [Axt,  Ax2,  •  •  •  ,  Ax„],  Ax  can  be  solved 
by  using  a  globally  convergent  modified  Newton’s  method 
[6],  and  the  parameter  vector  can  be  modified  by 

xKW  =  xM  +  Ax.  (9) 

Thus  the  optimum  solution  x*  of  S  is  further  approached 
by  xnew.  The  modified  Newton’s  iteration  can  be  written 
as 

Ax*+I  =  Ax*  -  X[e(Ax*)]"'g(Ax*)  (10) 
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where 


where 


0(Ax) 


dg(Ax) 

dAx 


=  p(p  -  \)BDBt  +  2 pV  +  2p{p  -  1  )G 

(11) 

where  V  -  diag  C2R ,  ■  •  •  ,  C*/?].  G  =  diag 

[C{DCtY,  C2DCtY,  •  ••  ,  CnDCrY],  and  X  can  be  found 
by  backtrack  linear  search  [6]  such  that 

S(Ax*  +  l)  <  S(Axk)  -  aX[g(Ax*)jr 

[Q(Ax*)]‘‘  •  g(Ax*)  (12) 

where  a  =  0.0001  [6].  Note  that  there  is  only  a  simple 
polynomial  calculation  rather  than  a  complicated  objec¬ 
tive  function  evaluation  required  for  the  linear  search  em¬ 
ployed  here,  which  relieves  the  computation  load  sub¬ 
stantially.  Moreover,  if  the  Hessian  matrix  Q (Ax)  = 
V2S(x)  is  nonpositive  definite,  we  can  convert  it  to  a  pos¬ 
itive  definite  one  by  substituting  Q(Ax)  by  [Q(Ax)  4- 
yl]  (y  >  0)  [6]  to  avoid  the  possible  convergence  at  the 
maximum  or  saddle  point. 

Once  convergence  of  Newton’s  iteration  is  reached, 
(i.e.,  xnew  could  be  found)  we  calculate  5new  =  S(xnew) 
and  compare  it  with  S'"1*  and  SM. 

If  Sn"“  <  Smax,  i.e.,  the  current  iteration  is  successful, 
we  update  the  super-polygon  by  replacing  xma*  with  x”'*; 
further  if  5new  <  SM,  we  update  the  super-polygon  by 
replacing  xmax  with  xM  and  substitute  xncw  for  xM.  After 


SD 


2n  +  l 

E 

k  =  1 
k  *  max 


-1 

5* 


Then  a  certain  reflection  point  of  xmax  with  respect  to  x  is 
chosen  as  the  new  point,  i.e., 

x**  =  X  +  0(x  -  x™x),  0e(O,  1],  (14) 

If  Snew  =  5(xnew)  <  SM  or  Snew  <  Smax,  the  same 
procedure  for  successful  iteration  mentioned  above  can  be 
resumed. 

If  SM*  >  Smax  still  holds,  we  can  apply  the  quadratic 
interpolation  scheme  through  the  three  known  points  xmax, 
xM,  and  x1**  to  find  a  new  point  Xne'1',  since  SM  <  Smax 
<  S'1'"  in  this  case.  If  the  value  of  the  objective  function 
at  X new  is  smaller  than  5max,  xnew  is  replaced  by  Xnrw.  In 
case  Sne *  =  S(Xnew)  is  still  greater  than  Smax,  which  rarely 
happens,  the  minimum  point  of  5  may  be  inside  the  super¬ 
polygon.  If  this  happens,  we  can  shrink  the  size  of  the 
super-polygon  to  locate  the  minimum  point  more  accu¬ 
rately  [2]. 

However,  based  on  our  experience,  the  odds  of  needing 
to  shrink  the  super-polygon  are  quite  rare.  As  a  result, 
POSM  needs  evaluations  of  objective  functions  no  more 
than  twice  for  any  failed  iteration  except  in  some  extreme 
cases. 

The  algorithm  as  stated  above  can  be  summed  up  in  the 
following  algorithmic  procedure: 

Algorithm  POSM2 

(1)  Initialize  a  super-polygon  from  the  2n  +  1  points 
picked  around  the  initial  point  xf. 


x?(  1.02  4-  O, 

k  -  i. 

x?  ^  0,  k  <  n 

(1-02  4-  O, 

k  =  i. 

x°  =  0,  k  <  n 

x?(0.98  -  $). 

(k  -  n)  =  t. 

xf  *  0,  n  <  k  <  2n 

(0.98  -  £), 

(k  -  n)  =  i, 

xf  =  0,  n  <  k  <  2n 

x° 

i  * 

k  *  /, 

(k  -  n)  *  i 

the  reference  point  xmax  is  determined  from  the  new  super¬ 
polygon,  the  next  optimization  iteration  starts  until  the 
convergence  condition  is  satisfied  Note  POSM  only  eval¬ 
uates  the  objective  function  once  in  the  successful  itera¬ 
tion. 

If  5new  >  5max,  the  current  iteration  fails.  The  remedy 
for  the  failure  case  is  a  crucial  point  in  this  algorithm.  If 
Snew  >  Smax  is  the  case,  we  apply  the  modified  simplex 
method  by  finding  the  weighted  center  x  of  the  super¬ 
polygon  with  2 n  vertices  (with  xmax  excluded).2 


x 


1 

SD 


2  n  +  I 

E 

k~  1 
k  *  max 


sk 


(13) 


JNoie  that  (he  weighied  center  of  the  super-polygon  is  chosen  as  new 
reference  point  in  MSM  employed  here,  instead  of  the  geometric  center  in 
conventional  SM  (2J. 


where  $  is  the  spreading  coefficient3  used  to  determine  the 
size  of  the  super- polygon. 

(2)  Evaluate  the  error  functions  e,  =  ( j  =  1 ,  2,  •  •  •  , 
m)  and  the  objective  functions  Sk. 

(3)  Select  the  point  with  minimum  error  xM  =  xmm 
among  2n  +  1  points  (SM  =  5(xM)}. 

(4)  Determine  the  point  with  maximum  error  xmax 
among  2n  +  1  points  {Smax  =  S(xmax)}. 

(5)  Calculate  pseudo  objective  function  S. 

(6)  Solve  Ax  using  a  globally  convergent  modified 
Newton  method  [6],  and  update  xnew  =  xM  4-  Ax. 

(7)  If  the  optimization  iteration  converges,  then  stop. 

(8)  If  Snew  <  SM,  then  update  the  super-polygon  by 
replacing  xmax  with  xM,  substitute  xnew  forxM,  and  go  to 
(4). 

’Selecting  1.02  and  0.98  instead  of  1.0  makes  points  jr  (t  =  I. 
2,  •  •  •  .  2n)  distribute  in  a  more  random  way. 
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(9)  If  Sne"  <  Smax,  then  update  the  super-polygon  by 
replacing  xmax  with  xnew  and  go  to  (4). 

(10)  (In  this  case,  Snew  >  Sma *  means  the  current  it¬ 
eration  fails.)  Find  xnew  by  the  modified  simplex  method 
(MSM)  and  calculate  Snew  =  S(x“w). 

(10.1)  If  Snew  <  SM,  then  update  the  super-polygon 
by  replacing  xmax  with  xM,  substitute  xnew  forxM,  and  go 
to  (4). 

(10.2)  If  Sncw  <  5max,  then  update  the  super-poly¬ 
gon  by  replacing  xmax  with  xnew,  and  go  to  (4) 

(10.3)  Calculate  a  minimum  point  Xnew  by  quadratic 
interpolation  through  the  three  points  (xnew,  Snew),  (xM, 
5"),  and  (*"“,  Smax),  then  replace  xnew  by  inew  and  cal¬ 
culate  Snew  =  5(xnew). 

(10.4)  If  Snew  <  Smax,  then  go  to  (8),  otherwise  go 
to  (11). 

(11)  Shrink  the  size  of  the  super-polygon,  i.e.,  calcu¬ 
late  (1] 


xk  =  x“  + 


(SM) 


h\Wp 


{Sk)'/P  +  (SM)'/p 


(x*  -  xM). 


(12)  Then  go  to  (2). 

Notes: 

(1)  In  step  (1),  users  only  need  to  provide  an  initial 
point  x°.  The  other  2 n  points  will  be  chosen  by  moving 
x  along  positive  and  negative  directions  on  each  coordi¬ 
nate  axis,  and  the  distance  is  determined  by  the  spreading 
coefficient.  Results  show  that  POSM  is  not  sensitive  to 
the  selection  of  the  initial  super-polygon. 

(2)  In  step  (7),  if  (1)  |SmiU  -  SM|  <  e,S"  +  €2  and 
(2)  |x,max  -  xf  |  <  e,xM  +  f,  (i  =  1,  2,  *  •  •  ,  n),  then 
the  iteration  stops.  Here,  e2  is  the  absolute  error  tolerance 
and  fi  is  the  relative  error  tolerance  which  are  given  by 
the  user  (e2  =  e2  =  10~12  by  default).  For  the  examp'es 
in  Section  III,  the  stopping  criterion  is  set  as  ||  g  |j  <  or 
|  S|  <  «|  in  order  to  have  fair  comparisons  with  the  results 
in  [1]  and  [7]. 

(3)  In  step  (10.3),  the  quadratic  interpolation  function 
is  applied  to  the  three  points  (xnew,  Snew),  (xM,  SM),  and 
(xmax,  5max)  and  yields 

»(x) .  *-  >  -y/y  -  - 

(xmax  -  xM)  (x™x  -  x"ew) 

,  ^  (*  -  -  xnn 

(xw  -  x™x)r(xM  -  xnew) 

+  5n.w  (x-X™X)r(x-x") 

(x^  -  xmax)r(xne'*  -  xM) 

The  minimum  point  of  Sl(x)  can  be  expressed  as 

new  =  (SmaxRC  4  S”**AB)xM  +  (SM/1C 

*  2  (Smax 


(4)  From  a  number  of  examples  in  the  following  sec¬ 
tion,  we  will  find  that  the  case  of  (11)  occurs  only  when 
xM  is  very  close  to  x*  (the  optimal  point).  In  general,  if 
this  happens,  unless  very  high  accuracy  is  needed,  the 
iteration  can  be  terminated  because  the  improvement  of 
1 5(x)  |  will  be  very  little  if  the  iteration  keeps  going. 

(5)  If  Chen’s  method  is  employed  to  calculate  Ax,  in 
which  a  linear  approximation  of  e;  and  a  second-order  ap¬ 
proximation  of  ${xM  +  Ax)  are  assumed,  and  the  steps 
(10)  and  (11)  are  adopted  for  processing  the  failed  itera¬ 
tion  instead  of  Chen’s  strategy  [1],  we  would  call  this 
modified  algorithm  the  first-order  POSM  (POSM1),  and 
the  algorithm  mentioned  above  the  second-order  POSM 
(POSM2)  (or  abbreviated  to  POSM). 

(6)  It  is  evident  that  the  POSM  needs  at  least  2 n  +  1 
initial  objective  function  calculations,  while  Chen  s 
method  needs  only  n  +  1  initial  objective  function  cal¬ 
culations.  In  order  to  circumvent  this  shortcoming  a 
variable  order  POSM  scheme  has  been  proposed  in  which 
the  first-order  POSM  is  applied  in  the  first  n  iterations 
without  discarding  the  replaced  points;  then  the  second- 
order  POSM  is  initialized  through  the  super-polygon 
formed  by  the  accumulated  2 n  +  1  points  automatically. 

III.  Numerical  Examples  and  Comparisons  with 
Other  Algorithms 

First,  we  will  present  the  computation  results  of  some 
numerical  examples  to  show  the  advantage  of  the  high 
efficiency  embodied  in  the  POSM  in  comparison  with 
other  existing  methods. 

Example  1:  Rosenbrock’s  Parabolic  Valley  problem 
(P  =  2): 

S  =  100! x,  -  x]f  +  (1  -  x,)2. 

The  objective  function  has  a  minimum  value  of  zero  at 
( 1,  1 ).  The  initial  point  is  set  at  (  -  1.2,  1 ).  The  results 
are  listed  in  Tables  I  and  II,  where  AS  is  the  change  of 
the  objective  function. 

As  seen  in  the  results,  the  POSM  is  superior  to  other 
well  known  methods  in  terms  of  a  minimum  number  of 
objective  function  evaluations.  Also,  it  is  an  obvious  ad¬ 
vantage  that  POSM  is  insensitive  to  the  spreading  coeffi¬ 
cient 

The  objective  function  and  the  POSM  optimization  pro¬ 
cess  are  illustrated  in  Figs.  1  and  2  by  showing  contours 
of  constant  values  of  objective  function  and  the  evaluated 
variable  loci  of  this  example  as  the  optimization  proceeds. 
It  is  clearly  found  that  the  variable  loci  move  downward 


+  Snew/4B)xmax  +  ( Sm"BC  +  SM/fC)xne'v 
+  S,w  +  Snew) 


where 

A  =  (xmax  -  x'V(xmax  - 

B  =  {xM  -  xmax)r(xM  -  xm) 
C  =  (xm  -  xmax)r(xnew  -  xM) 


towards  the  local  bottom  of  the  objective  function  valley 
at  the  first  step,  then  continuously  move  along  the  mini¬ 
mum  objective  function  valley  toward  the  optimal  point 
with  fast  monotonic  decrease  in  objective  function  value. 

Example  2:  Colville’s  four-dimensional  banana  prob- 
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Fig.  1.  The  objective  function  contours  of  Example  1  showing  the  fast 
monotonic  decrease  in  objective  function  value  at  consecutive  optimi¬ 
zation  steps  by  POSM;  OFV  =  5,  OFV  =  10.  •  •  • :  OFV  =  15, 
OFV  =  20. 


Fig.  2.  A  three-dimensional  plot  of  variable  loci  illustrating  the  optimi¬ 
zation  path  of  Rosenbrock's  parabolic  valley  problem  by  POSM,  and  * 
symbols  indicating  the  initial  super-polygon. 


TABLE  I 


Method 

Stopping  Criterion 

Ns  (Nsg)  t 

Fletcher-Powell 

S 

< 

icru 

120  -  380*[ll 

Jacobsoo-Olcsman 

A 

s 

<  1(T* 

(39)  (I] 

Powell 

S 

< 

1(T* 

70  [1] 

Peckham 

S 

< 

4xlO~10 

12(1] 

Hoofce-/eevs 

s 

< 

2  xio-7 

353  (71 

Nelder-Mead 

s 

< 

10’7 

300  [7J 

Powell 

s 

< 

10'M 

150  -  1500*[7] 

Newton 

s 

< 

to-'3 

98(7] 

D-S-C 

s 

< 

ur7 

187  [7] 

Simplex 

s 

< 

ter20 

133-161* 

Cbeo  [1] 

s 

< 

lcr20 

7-135* 

POSM 

s 

< 

10-» 

15-17* 

•Spreading  coefficients  ({ )  0. 1 . 0.4,  1.0,  2.0  are  used  for  these  results. 
$Ns  is  the  total  number  of  objective  function  evaluations;  Nsg  is  the 
total  number  of  objective  function  calculations  and  its  associated  deriva¬ 
tives. 


TABLE  II 


Algorithm 

4-0.1 

0.4 

!  1.0 

2.0 

Simplex 

161 

137 

133 

140 

Chen 

68 

135 

71 

7 

POSM 

16 

15 

17 

16 

lem  ( p  =  2) : 

5  =  100(x?  -  x2f  +  (1  -  x,)2  +  90(^3  -  xj 
+  (1  -  x2f  +  10.1  [(x2  -  l)2  +  (x4  -  l)2] 

+  19.8(x2  -  l)(x4  -  1). 

Note  that  this  objective  function  is  not  an  explicit  sum 
of pth  error  terms  in  real  numbers,  but  it  can  be  rearranged 
in  a  least  pth  approximation  form  with  p  =  2  like  [1] 

5  =  2  ej 

i- i 

where  the  error  terms  are  defined  as 
e,  =  10(x?  -  x2) 

*2  =  (1  ~  *:) 
e3  =  V90(x]  -  x4 ) 

*4=0  -  X2) 

e5  =  V02(x2  -  1) 

*6  =  V0^(X4  -  1) 
e-j  =  V9.9(x2  +  x4  —  2). 

The  objective  function  has  a  minimum  value  of  zero  at 
( 1 ,  1 ,  1 ,  1 )  and  a  stationary  point  at  (  -0.9679,  0.947 1 , 
—0.9695  ,  0.9512).  The  evaluation  results  of  POSM  with 
various  initial  points  are  summarized  in  Table  III  for  com¬ 
parison  with  other  methods. 

It  is  interesting  to  point  out  that  the  POSM2  algorithm 
has  a  convergence  rate  faster  than  that  of  any  other  meth¬ 
ods  even  when  the  initial  point  happens  to  be  the  nonop- 
timal  stationary  point.  In  addition,  the  results  show  that 
the  location  of  the  initial  point  does  not  affect  the  conver¬ 
gence  rate  much  in  POSM2. 

Example  3:  Powell's  Function  problem  (p  =  2): 

S  -  (x,  +  10x2)‘  +  5(x3  -  x4)‘ 

+  (x2  -  2x3)J  +  10(x,  -  x4)4. 

The  objective  function  has  a  minimum  value  of  zero  at 
(0,  0,  0,  0)  where  the  Hessian  matrix  is  singular.  Table 
IV  shows  the  results  calculated  by  POSM  in  contrast  with 
those  obtained  by  other  methods. 

Even  though  the  POSM  does  not  employ  an  exact  sec¬ 
ond-order  formulation  by  neglecting  the  mixed  terms  for 
reducing  the  complexity  of  objective  function  evaluation 
in  the  error  function  approximation,  the  performance  of 
the  POSM  is  still  superior  to  other  methods,  since  the  rc 
suits  from  the  preceeding  examples,  wnere  there  are  mixed 
terms  involved  only  in  the  error  functions  e2  and  e4  of 
Example  3,  show  that  the  inclusions  of  the  mixed  terms 
in  the  error  function  has  only  a  minor  influence  on  the 
convergence  rate.  Thus  we  trade  the  accuracy  in  the  error 
function  approximation,  where  higher  accuracy  in  ap¬ 
proximation  may  lead  to  faster  convergence,  for  a  less 
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TABLE 

min 

Method 

Jacobson 

-Ofcunao 

Hgjsio-* 

Fletcher 

-Powell  | 
Uglisio-1 

Simplex 

Chen 

POSM 

Initial  Point 

ASS  KT* 

M  s  10-* 

ng!i  s  to-* 

(-3.-l.-3.-l) 

135 

161 

694 

26 

25 

(-3.1.-3.I) 

154 

w 

549 

31 

26 

(-uu.-i.2j) 

194 

64* 

573 

107 

30 

sodoouy  point 

fail 

fail 

3*2 

231 

41 

TABLE  IV  [1] 


Method 

Fletcher 

-Powell 

llgusur* 

Jacobson 

-Ckjsman 

llgpio-* 

Rank-one 

Project- 

Huang  ! 

Chen 

POSM 

Initial 

Point 

Ugli^xo-6 

IgllSHr4 

Ilglisio-4 

IlgllSIO-* 

"glisicr4 

3.-1.0.1 

*0 

64 

28  , 

21 

10.10.10,-10  j 

378-925 

286-660 

172-663 

28  j 

43 

computation  load  of  objective  function  evaluation  in  order 
to  maximize  the  overall  performance. 

IV.  Circuit  Design  Examples 

Due  to  its  inherent  high  efficiency  advantage,  POSM  is 
highly  suitable  for  applications  in  the  optimization  of 
complicated  integrated  circuit  designs.  As  a  rule,  the  more 
complex  the  IC  is  the  more  efficient  it  appears. 

Example  4:  AC  optimization  of  a  second-order  active 
low-pass  filter: 

The  circuit  schematic  is  shown  in  Fig.  3  and  the  input 
circuit  description  codes  are  listed  in  (1)  of  Appendix  2. 

Fig  4  illustrates  the  frequency  response  of  the  circuit 
at  different  optimization  stages.  As  seen,  a  significant  dis¬ 
crepancy  exists  between  the  corresponding  frequency  re¬ 
sponse  (#2)  for  the  initial  parameter  vector  and  the  de¬ 
sired  frequency  response  (#1).  In  a  few  iterations,  the 
intermediate  optimized  response  (#3)  has  been  very  close 
to  the  desired  response.  Finally,  a  fairly  good  agreement 
between  the  desired  response  and  the  optimized  response 
(#4)  is  observed  with  reasonable  optimized  parameter 
values.  The  quantitative  details  that  include  iteration  sta¬ 
tistics,  objective  function  values,  and  parameter  vectors 
are  listed  in  Table  V  and  compare  with  the  results  ob¬ 
tained  by  Chen’s  method  [1],  (run  on  IBM  4381,  the  same 
hereafter)  where  the  performances  are  evaluated  in  terms 
of  the  following: 

A/,  the  number  of  iterations, 

Nc  the  number  of  objective  function  evaluations, 

,V„  the  number  of  additional  objective  function  evalu¬ 
ations  while  iteration  fails,  ?rd 
IVOF,  FVOF,  IPV,  FPV  stand  for  initial  value  and  fi¬ 
nal  value  of  objective  function,  initial  and  final  pa¬ 
rameter  vector,  respectively, 
p  is  the  order  of  least  pth  approximation. 

From  Table  IV,  it  is  found  that  Chen’s  method  is  sen¬ 
sitive  to  p  (large  variation  on  the  number  of  iterations 


Fig.  3.  The  circuit  schematic  of  an  active  low-pass  filter. 


Fig.  4.  The  simulated  frequency  responses  of  the  active  low-pass  filter 
shown  in  Fig.  3.  #1  the  desired  ac  response;  #2  #4  the  initial  calculated 
frequency  response  and  the  improved  frequency  responses  obtained  at 
successive  optimization  stages. 

with  respect  top),  while  POSM  is  not,  and  POSM  is  more 
effective  than  Chen’s  no  matter  what  the  p  value  is. 

Moreover,  the  final  row  in  Table  V  shows  if  linear 
search  (as  in  Chen’s  method)  is  adopted  instead  of  MSM 
in  step  (10)  in  the  algorithm  POSM2,  then  the  number  of 
iteration  and  failed  iteration  increases.  It  is  apparent  that 
the  low  efficiency  inherent  in  the  linear  search  strategy 
causes  the  degradation  in  the  performance  of  Chen’s 
method. 

Example  5:  Parameter  extraction  of  MOS  device: 
Because  of  technological  limitations,  the  final  values  of 
the  devices’  parameters  will  not  match  the  design  values 
exactly  especially  for  small  size  devices.  In  order  to  pre- 
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TABLE  V 


Method 

Hr 

Nc 

Na 

P 

IVOF 

FVOF 

tPV 

FPV 

CPU  I:n.c 

Chen 

43 

211 

165 

2 

3.120 

xlO3 

429 

Cl-lOU 

C2-IU 

CMJ82U 

C2-0.315U 

83.37  S 

POSM 

20 

25 

0 

2 

3.120 

xlO3 

0  0987 

CU10U 

C2»tU 

CU3.02U 

C2=O.097U 

12.37  S 

Chen 

16 

23 

4 

4 

1.736 

xlO* 

0.082 

C1-10U 

C2-1U 

CU3.17U 

C2^).095U 

11  30  S 

POSM 

. 

14 

19 

0 

4 

1.736 

xlO* 

0.071 

CI-10U 

C2-IU 

C1-3.16U 

C2-0.092U 

9  870  S 

Chen 

34 

77 

40 

6 

4.i6; 

xIO* 

23.9a 

Cl-.GU 

C2»lO 

Ci-3  48U 

C2-43  083U 

31.34  S 

POSM 

21 

28 

2 

6 

4.169 

xlO’ 

0.39 

C1-10U 

C2=1U 

CU3.25U 

C2=0.094U 

13.30  S 

POSM* 

24 

33 

4 

6 

4.169 

xlO9 

041 

CU10U 

C2*IU 

Cl  *3. 29  U 

C2  =<0.09111 

*4.05  5 

•For  the  handling  of  iteration  failure  here.  Chen's  method  is  used  rather 
than  the  Modified  Simplex  Method  as  in  step  (10)  of  POSM2. 


pare  a  more  accurate  empirical  device  model,  the  devices’ 
effective  parameters  (e.g.,  channel  length)  must  be  ex¬ 
tracted  from  measured  data.  To  this  end,  the  output  char¬ 
acteristics  of  device  are  measured  first  and  taken  as  the 
given  requirement;  then  the  dc  optimization  is  performed. 

The  circuit  for  parameter  extraction  of  a  MOS  device 
is  depicted  in  Fig.  5,  and  its  input  codes  are  listed  in  (2) 
of  Appendix  2.  The  results  as  compared  with  those  ob¬ 
tained  by  Chen’s  method,  shown  in  Table  VI,  show  that 
if  the  initial  value  L  =  10(7  is  selected,  the  number  of 
iterations  by  Chen’s  method  increases,  while  the  number 
of  iterations  in  POSM  does  not  change  whether  L  =  9 U 
or  L  =  10 U.  The  results  of  POSM  are  not  sensitive  to  the 
initial  value  of  parameters. 

Example  6:  Transient  characteristic  optimization  of  a 
circuit  simulating  a  low  power  position-follow  system: 

In  this  example,  the  objective  is  to  make  the  rise  and 
fall  times  of  the  TR  response  shorter  and  its  overshoot 
smaller.  The  results  compared  with  those  obtained  by 
Chen’s  method  are  shown  in  Table  VII,  and  the  transient 
characteristics  in  Figs.  7-10.  The  circuit  schematic  is 
shown  in  Fig.  6. 

As  shown  in  Table  VII,  POSM  comes  up  with  a  better 
performance  than  Chen’s  method,  which  has  very  often 
been  involved  with  linear  searching  and  simplex  shrink¬ 
ing,  and  the  number  of  additional  objective  function  eval¬ 
uations  is  enormous.  Therefore,  POSM  is  preferable.  Ini¬ 
tially  in  this  example,  the  circuit  was  optimized  using 
POSM  with  weight  as  1  (default  value).  The  output  wave¬ 
form  of  this  optimization  as  shown  in  Fig.  8  did  not  meet 
the  design  goal;  therefore,  the  optimization  was  further 
carried  out  interactively  by  altering  the  weights  (i.e.,  re¬ 
ducing  the  weight  at  the  sampling  point  where  the  devia¬ 
tion  is  small  and  increasing  it  at  the  point  where  the  de¬ 
viation  from  the  desired  value  is  large).  For  example,  the 
weights  of  2nd-6th  and  1 3th-  17th  points  are  raised  to  10. 
The  optimization  process  is  stopped  once  the  output  char¬ 
acteristic  comes  closest  to  the  desired  response.  The  re¬ 
sults  are  shown  in  Fig.  9.  The  initial  calculated  voltage 
waveform  of  w (21,  0)  in  the  circuit  as  shown  in  Fig.  o  is 


(7) 


"VIDI>  C7  VU>‘ 

.Jq 


© 


Fig.  5.  The  MOS  circuit  used  for  effective  device  parameter  extraction. 


Fig.  6.  The  equivalent  circuit  of  a  low  power  position-follow  system. 


TABLE  VI 


Method  j 

N, 

[  Nc 

N«  | 

IVOF  | 

l  FVOF 

IPV 

FPV  , 

CPU  nme 

Cben  | 

8 

Luj 

1  |  1.5287 

4JXIQ-’ 

L=9U 

L-6.00U 

5  86  S 

POSM  j 

■  4  1 

l  °  I 

|  1.5787  J 

4.3x10-’ 

L-9U  j 

L-6.00U  , 

|  4  76S 

Chen  | 

9  | 

12 

1  |  1.5787  | 

4.3x10-’ 

L-10U 

U6  00U 

|  7.06  S 

POSM 

4 

7 

0  i 

1.5287 

4.3x10-’ 

L»10U  j 

L-6.00U  | 

|  4.24  S 

TABLE  VII 


Method  j 

Nt 

Nc 

Na  j 

!  IVOF 

FVOF  | 

IPV  | 

FPV 

CPU  tune 

;  i 

R~-  -00k 

R2-4I3K  j 

Chen 

is 

76 

48 

0.2548 

0.2398 

R3*l00k 

R3*92K  | 

1836  S 

R3*50k 

R4-507K  j 

i 

R 4= 100k 

R5*126K 

POSM 

R7»400k 

R7*4l 1  IK  1 

with 

18 

27 

0 

02548 

0.067 

R  3*1 00k 

R3-97.7K  | 

575  S 

uniform 

R4*50k 

R4.313K 

weight  1  | 

1 

R5*l00k  ( 

R5-996K  j 

POSM 

R2»400k 

R2-679K 

with 

17 

26 

0 

7.918 

0.2663 

R3*100k  1 

R3-I28.6K 

389  S 

vanable  : 

I 

R4»50k 

R4-16.2K  ! 

weight!  | 

R5*l00k 

R5-86  6K  | 

illustrated  in  Fig.  7.  The  waveform  obtained  by  Chen’s 
method  is  shown  in  Fig.  10.  It  is  evident  that  the  over¬ 
shoot  as  well  as  the  rise  and  fall  times  are  reduced  re¬ 
markably  after  POSM  optimization  application. 

Example  7:  Optimization  of  a  high-order  bandpass  fil¬ 
ter  design 

As  shown  in  Fig.  11,  the  circuit  consists  of  an  active 
7-order  low-pass  filter  and  an  active  5-order  high-pass  fil¬ 
ter.  The  specification  of  the  design  requires  that  the  atten¬ 
uation  in  stop-band  should  be  greater  than  45  dB.  There 
are  22  design  parameters  assigned.  The  input  circuit  de- 
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OUTPUT  *  -OSV  00V  05V  1,0V  1.5V 

TM6=S 
ooooe.oo 

I  0006-01 
20006-01 
30006-01 
4  0006-01 
50006-01 
6  0006-01 
70006-01 
8  0006-01 
90006-01 
1  0006.00 
1  1006*00 
1  2006.00 
1  3006.00 
1  4006.00 
1,5006.00 
1  6006.00 
1  7006.00 
1  8006.00 

1  9006.00 

2  0006,00 


OUTPUT  „  -OSV  0  0V  0.5V  10V  15V 

Tiwe-s 
0  0006.00 
1  0006-01 
2  0006-01 
30006-01 

4  0006-01 

5  0006-01 
6.0006-01 

7  0006-01 

8  0006-01 
9.0006-01 
1,0006.00 
1. 1006.00 
1,2006.00 
1  3006.00 
1  4006.00 
1.5006*00 
1  6006.00 
1.7006.00 
1  0006.00 

1  9006.00 

2  0006.00 


Fig.  7.  The  initial  voltage  transient  response  at  o(0,  21)  of  the  circuit 
shown  in  Fig.  6. 


Fig.  9.  The  optimized  transient  response  at  t>(0,  21 )  of  the  circuit  shown 
in  Fig.  6  using  the  POSM  with  variable  weights. 


OUTPUT  =  -0.5V  0.0V  0.5V  1.0V  1.5V 

T1M6.S 
0  0006,00 
1.0006-01 
2  0006-01 

3  0006-01 

4  0006-01 

5  0006-01 

6  0006-01 
7.0006-01 
8.0006  01 
9  0006-01 
1  0006.00 
1  1006.00 
;  2CCP • 00 
1  3006,00 
1  4006,00 
I  5006,00 
1  6006,00 
1  7006.00 
1  8006,00 
1.9006,00 
20006,00 


OUTPUT  = 
TW6-S 
0.0006.00 
1.0006-01 
2.0006-01 
3.0006-01 
4.0006-01 
5.0006-01 
6.0006-01 
70006-01 
8.0006-01 
90006-01 
1.0006*00 
1.1006,00 
1.2006,00 
1.3006,00 
1  4006.00 
1.5006,00 
16006,00 
1  7006,00 
1.8006,00 
1  9006.00 
20006.00 


-0.5V  0.0V  0.5V  1  OV  1  5V 


Fig.  8.  The  optimized  transient  response  at  v(0,  21 )  of  the  circuit  shown 
in  Fig.  6  using  the  POSM  with  uniform  weights  l's. 


Fig.  10.  The  optimized  transient  characteristic  calculated  by  Chen's 

Method. 


scription  codes  are  listed  in  (4)  of  Appendix  2.  The  opti¬ 
mization  results  are  shown  in  Table  VIII. 

IPV for  both  are: 

R4  =  100K  R15  =  100K  R28  =  100K 

R5  =  100K  R19  =  100K  R29  =  91K 

R6  =  80K  R20  =  100K  R30  =  50K 

R7  =  10K  R21  =  50K  R34  =  100K 

R8  =  50K  R22  =  100K  R35  =  81K 

R12  =  100K  R26  =  100K  R36  =  10K 

R13  =  100K  R27  =  100K  R37  =  50K 

R14  =  50K 


FPV for  Chen 's  method  are: 

R4  =  157. 8K  R15  =  -36. 3K  R28  =  176. 6K 

R5  =  144. 8K  R19  =  36.52K  R29  =  -47.0K 

R6  =  4 1.2  IK  R20  =  64.57K  R30  =  25.24K 

R7  =  8.35K  R21  =  56.05K  R34  =  287.5K 

R8  =  46.45K  R22  =  34.57K  R35  =  82K 

R12  =  20.6K  R26  =  194. 8K  R36  =  -9.09K 

R13  =  98.32K  R27  =  319.2K  R37  =  159.7K 

R14  =  59.03K 


FPV  for  the  POSM  are: 

R4  =  1 18. 7K  R15  =  30.51K  R28  =  104.9K 

R5  =  100.6K  R19  =  126.0K  R29  =  84.59K 

R6  =  92. 1  IK  R20  =  124.  IK  R30  =  31.79K 

R7  =  7.71K  R21  =  8.64K  R34  =  98.98K 

R8  =  49.97K  R22  =  90.74K  R35  =  82.7  IK 

R12  =  24.34K  R26  =  90.83K  R36  =  24.99K 

R13  =  94.65K  R27  =  97.86K  R37  =  66.44K 

R14  =  1000K 

Notice  that  the  results  shown  in  Table  VIII  by  the 
POSM  were  summarized  from  a  convergent  computer  run, 
whereas  those  by  the  Chen's  method  were  extracted  from 
intermediate  optimization  results,  which  are  not  accept¬ 
able  due  to  unrealistic  resistances  because  it  fails  to  con¬ 
verge  after  a  long-time  computer  run. 

This  fact  shows  that  for  large  size  optimization  prob¬ 
lems,  the  POSM  basically  runs  in  time  proportional  to  the 
dimensions  of  the  parameter  vector  and  always  obtains 
convergent  results,  while,  on  the  contrary,  Chen's  method 
often  fails  to  reach  a  reasonable  solution  because  it  easily 
becomes  trapped  in  the  linear  searching  and  simplex 
shrinking  loop. 
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TABLE  VIII 


Method  ] 

Ni 

Nc  | 

Na  j 

IVOF 

FVOF 

CPU  Tune 

POSM  ' 

56 

101 

0 

16350 

4562 

560  S 

Chen  | 

* 

131 

_ 1 

707 

553 

16350 

1 

5655 

2645  S 

1 _ - 

Example  8:  Constrained  optimization  on  the  transient 
response  of  an  inverter  chain  [8]: 

In  this  example,  the  signal  delay  is  subjected  to  min¬ 
imization  under  the  constraint  that  transient  power  is  less 
than  2  mW.  The  channel  lengths  of  PMOS’s  and  widths 
of  NMOS’s  are  designated  as  adjustable  parameters  with 
a  total  number  of  10.  The  circuit  schematic  is  shown  in 
Fig.  12.  The  results  are  shown  in  Table  IX  and  the  initial 
and  optimized  output  responses,  where  the  signal  delay  is 
reduced  from  td  =  148  ns  to  td  =  80  ns,  are  plotted  in 
Figs.  13-14,  respectively. 

V.  Conclusions 

We  have  presented  an  unconstrained  optimization  al¬ 
gorithm  suitable  for  integrated  circuit  design  applications. 
The  algorithm  has  been  implemented  in  the  general-pur¬ 
pose  analysis  and  design  program  for  integrated  circuit 
application  ADIC.  To  show  the  effectiveness  and  effi¬ 
ciency  of  this  algorithm,  we  have  evaluated  it  by  a  num¬ 
ber  of  typical  numerical  examples  which  are  commonly 
examined  by  other  existing  optimization  methods  as  well 
as  some  practical  IC  designs.  As  a  result,  the  POSM  has 
proved  to  be  more  effective  and  efficient  than  other  meth¬ 
ods  in  terms  of  the  number  of  iterations,  objective  func¬ 
tion  evaluations,  the  chance  of  needing  super-polygon 
shrinkage,  and  CPU  time.  The  high  efficiency  is  attrib¬ 
uted  to  1)  POSM  adopts  high  order  approximations  for 
error  function,  2)  the  modified  simplex  method  is  used  for 
handling  the  failed  iteration  during  optimization,  and  3)  a 
globally  convergent  modified  Newton  method  is  used  for 
obtaining  the  minimum  point  of  S.  By  this  sort  of  com¬ 
bination,  POSM  reduced  bcth  the  number  of  objective 
function  calculations  and  the  chance  of  shrinking  poly- 


Fig.  12.  The  circuit  schematic  of  an  inverter  chain. 


OUTPUT  -  0  000E+00V  2.D00E+00V  4.0QOE+OOV  6  OOOE+OOV 

TIMERS  - 
0.00QE+00S 
1.000E-08S 
2.0006-08$ 

3  00QE-08S 
4.000E-08S 
5.000E-08S 

6  000E-08S 

7  000E-08S 

8  0006-08S 

9  0006-08S 
1.0006-07$ 

1.1006-07S 
1.2006-07S 
1.300E-07S 
1.400E-07S 
1  500E-07S 
1  600E-07S 
1.700E-07S 
1.80C6-07S 
1  90CE-07S 
2000E-07S 


Fig.  13.  The  initial  output  response  of  an  inverter  chain  (/rf  =  148  ns). 


OUTPUT  =  0  00CE*00V  2.000E+00V  4  0C0E*00V  6  OOOE+OOV 


Fig.  14.  The  optimized  output  response  of  an  inverter  chain  (rrf  *  80  ns). 
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TABLE  XI 


Method 

rt 

"Re 

rJ  a 

IVOF 

FVOF 

1PV 

FPV 

CPU  ume 

POSM 

34 

13 

0 

0.714x10* 

0.159x10s 

M1.L-20U 

L-1129U 

157.63  S 

M2.W-16U 

W.23.92U 

M3J_*20U 

L-5.360U 

M4.W-I6U 

W-16.63U 

M5.L-20U 

U17.58U 

M6.W-16U 

W-11.60U 

M7JL-20U 

L-M.9IU 

MS.W-16U 

W-I5.99U 

M9JL-20U 

L-19  54U 

M10.W-16U 

W-I5.36U 

gon,  and  led  it  insensitive  to  the  initial  polygon  and  the  p 
value.  Experiments  in  IC  design  applications  also  show 
that  the  POSM  performs  equally  well  in  terms  of  accuracy 
and  effectiveness  for  large  size  problems,  especially  for 
complex  circuits  with  very  complicated  objective  func¬ 
tions,  and  very  fast  in  computer  run-time.  Extensions  to 
constrained  optimization  have  also  been  proposed  and 
tested  for  emerging  high  density  VLSI  circuit  designs,  on 
which  the  design  constraints  imposed  are  more  stringent. 


EPS  is  the  accuracy  of  optimization  (default  0.001 ),  P  is 
the  order  of  the  least  pih  approximations  (default  2). 
CONS  specifies  the  constraints  where  C  can  be  PW 
(power),  TR  (rise^time),  TF  (fall  time)  and/or  TD  (delay 
time),  etc.,  and  C,  is  upper  limit  value  of  C,.  The  objec¬ 
tive  function  of  the  circuit  optimization  problem  can  be 
written  as 

k  nu 

S(x)  =  S  (^(«,(x,  !,)  -  (Al) 

where  u,(x,  i;)  is  the  calculated  value  at y'th  sampling  point 
of  Vj. 

For  the  constrained  optimization  problems,  the  La¬ 
grange  multiplier  method  (PHR  Algorithm  [10])  can  be 
applied.  In  general,  the  constrained  optimization  problem 
can  be  decomposed  into  a  series  of  unconstrained  opti¬ 
mization  problems4  whose  objective  function  is5 

/ 

Sc(x,  n,  X)  =  S(x)  +  E  X*  max  (l,(x),  -  X*//} 


Appendix  1 

Input  Statements  Used  in  ADIC-2  for 
Optimization  Purpose 

The  POSM  has  been  implemented  in  a  general-purpose 
IC  Analysis  and  Design  Program  ADIC-2. C.  The  input 
language  of  ADIC-2. C  is  compatible  with  that  of  SPICE- 
2G.5  with  some  additional  statements  included,  for  ex¬ 
ample,  parameter  statement,  optimization  statement, 
macro-model  statement  [8],  •  •  •  etc,  making  the  specifi¬ 
cation  of  optimization  problem  more  precise  and  com¬ 
plete. 


+  0.5/x‘  Z  [max  {l,(x),  -X*//x*}]' 

(A2) 

where  L,(x)  =  C,  -  C,  or  C,  -  C,,  p  is  the  penalty  func¬ 
tion  factor,  and  X,  is  the  Lagrange  multiplier  for  ith  con¬ 
straint.  It  should  be  noted  that  the  lower  and  upper  limit 
values  of  the  designable  parameters  are  treated  like  con¬ 
straints.  The  details  of  the  usage  are  described  in  Appen¬ 
dix  2. 


.PARA  [subcircuit-name,]  elementname  [(lower  limit, upper  limit)]  •  •  • 

+  [subcircuit-name,]  model-name, parameter-name  [(lower  limit, upper  limit)] 

.MIN  type  OBJT(Vl,  V2,  •  •  •  ,  Vk) 


isT 

11 

+ 

y.il.H'..] 

M.w'ul  •• 

+  w  =  w2 

VixlW. 2.) 

V2  2[,Wn)  •• 

+ 

+  w  =  wk 

Vk  ,[,H\,] 

Vkz[,Wk2]  ■■ 

■* 

+ 

% 

1! 

Srr 

t" 

*s 

II 

EPS  =  es. 

P  =  Mp, 

CONS{C\  (C, )  •  • 

■  •  (C/(C,)) 

where  .PARA  statement  is  used  to  specify  the  adjustable 
parameters  in  the  optimization  procedure,  which  could  be 
any  parameter  of  the  element  in  the  main  circuit  or  sub- 
circuits.  .MIN  statement  is  used  to  specify  the  optimiza¬ 
tion  objectives:  the  “type”  could  be  dc,  ac,  or  TRAN; 
the  optimization  objects  (OBJT)  indicates  the  nodal  or 
branch  variable  V,  to  be  fitted,  where  K,  could  be  branch 
currents,  nodal  voltages,  power  of  branches,  rise,  fall  or 
delay  time  of  signal  at  some  nodes,  and  W(,  Wxi  and  VXi 
are  the  overall  weight,  the  weight  and  desired  value  at yth 
sampling  point  of  variable  V„  respectively.  H^’s  are  Is 
by  default.  XI  is  the  spreading  coefficient  (default  0.03). 
LT  is  the  limit  of  the  number  of  iterations  (default  100). 


Appendix  2 

The  SPICE-Compatible  Circuit  Description  Codes 
for  the  Optimization  Examples  of  Integrated 
Circuit  Design 

(1)  A  Second-Order  Active  Filter 
A  second-order  active  filter 

.WIDTH  IN  =  80  OUT  =  80 
R1  I  2  283 
R2  2  3  729 
R3  0  4  10K 

*The  details  about  how  to  handle  the  constraints  will  be  presented  in 
another  paper. 

5C,  is  the  upper  limit  value  of  C,,  and  C,  is  the  lower-limit  value  of  C,. 
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R4  4  5  10K 
Cl  2  5  10U 
C2  3  0  1U 

XI  4  3  7  6  5  MACRO  OP-UA741 

VCC  6  0  DC  15 

VEE  7  0  DC  -15 

V1N  I  0  AC  1 

.AC  DEC  5  100  10K 

.PRINT  AC  VM(5)  VDB(5)  VDB(3) 

.PLOT  AC  VDB(5) 

.MIN  AC  OBJT(VDB(5))  W  =  1  6.054  6.054  6.054  6.054  6.054  6.054, 
0.1 

+  -0.3955  -7.960  -16.00  -24,05  -32.07  XI  =  0.5  EPS=0.0001 
.PARA  C1(100P.20U)  C2(100P,10U) 

.END 

In  this  example,  the  adjustable  parameters  are  Cl,  C2 
(.PARA  statement).  .MIN  statement  designates  ac  opti¬ 
mization  ;  the  output  to  be  optimized  is  V{  5 ) ;  there  are  1 1 
sampling  points6  that  are  determined  by  the  ac  statement, 
and  the  weight  of  every  sample  point  is  default  value  1 
except  the  sixth  point  whose  weight  is  0. 1 . 

The  macro-model  [11]  of  /xA741  is  included  for  mod¬ 
eling  the  operational  amplifier  in  this  circuit.  XI  ... 
MACRO  OP-UA471  statement7  is  used  for  calling  in  the 
macro-model  of  the  /xA741  from  the  macro-model  library 
and  its  parameters  can  be  modified  during  the  optimiza¬ 
tion  if  they  are  specified  in  the  .PARA  statement. 


C2  5  6  1U 
C3  7  nu 
C4  9  10  1U 
C5  17  19  1U 

C6  18  0  0.1U 

Cl  20  21  1U 

VIN  1  0 

+  PULSE(0  1  0  0.10  0.10  1.0  2.2) 

El  2  0 

+  POLY(2)  1  0  21  0  0  195  195 
E2  0  6  4  0  1 E 1 0 
E3  0  10  8  0  IE10 
E4  16  0  11  0  4.75 

E5  19  0  18  0  I 

E6  0  21  20  0  1 E 10 
R6  10  11  100K 

R7  11  12  10 

R8  11  14  10 

R9  11  0  100K 

V2  0  13  DC  6 
V3  15  0  DC  6 
D1  13  12  DWY 

D2  14  15  DWY 
.MODEL  DWY  D 
.TRAN  0.10  2.0  UIC 
.'LOT  TRAN  V(0,21) 

.PRINT  TRAN  V(2)  V(4)  V(6)  V(8) 

+  V(ll)  V(19)  V(21) 

.MIN  TRAN  O8JT(V(0,21))  W=1  0  1  1  1 
+  11111111000000000 
+  Xl  =  0.5  EPS  =0.001 
.PARA  R2  R3  R4  R5 
.END 


(2)  MOS  Output  Characteristics 


(4)  A  High-Order  Filter 


MOS  OUTPUT  CHARACTERISTICS 
.WIDTH  IN  =80  OUT  =  80 
VDS  3  0 
VGS  2  0 

Xll  120  SUBTEST 1 

.SUBCKT  SUBTEST  1  1  2  4 

Ml  1  2  4  4  MODI  L=  10U  W  10U 

+  AD  =  10P  AS=  10P 

.MODEL  MODI  NMOS  (VTO=-1.5 

+  UO  =  550  KP  =  3.0E-5  NSUB=1.E15 

+  LAMBDA  =0.02  PHI  =0.65) 

.ENDS 
V1D1  3  1 

.DC  VDS  0  8  .5  VGS  0  4  1 
.PRINT  DC  I(VIDI)  V(2) 

.PLOT  DC  I(V1DI) 

.MIN  DC  0BJT(I(V1DI) 

+  W=1000  0.0  3.16E-5  ...  (total  85 
+  ideal  values  at  sampling  points)  ... 

+  XI=0.5  EPS  =  0.0001 
.PARA  Xll ,Ml,L(2E-6,20E-6) 

.END 


(3)  Transient  Characteristics  Optimization  Example 

TRANSIENT  ANALYSIS  EXAMPLE 

R1  1  0  IK 

R2  2  3  400K 

R3  4  5  I0OK 

R4  6  7  50K 

R5  8  9  I00K 

R10  16  17  1.1571MEG 

RU  17  18  2.8692K 

RI2  19  20  1000K 

R40  4  0  1E8 

R80  8  0  1E8 

Cl  3  4  1U 


‘Here,  the  sample  points  are  the  output  points  of  response  K(5). 

7In  ADIC-2.C,  the  subcircuit  calling  statement  can  have  real  parame¬ 
ters,  i.e..  "XYYYYYYY  N1<N2  N3...)SUBNAM  [REAL-PARAME¬ 
TER]’. 


A  HIGH-ORDER  FILTER 

Cl  1  2  30U 

RI  3  2  220K 

R2  2  0  91K 

R3  4  40  15K 

Q1  3  2  4  QMOD1 

•••FIRST**** 

R4  4  5  100K 
R5  5  6  I0OK 
R6  6  9  80K 
R7  7  8  I OK 
R8  8  0  50K 
R9  II  0  10K 
RIO  7  11  51 
RI  1  3  10  10K 
C2  5  0  280P 
C3  5  8  230P 
C4  6  7  250P 
C5  8  9  1300P 
C6  8  0  HOP 
Q2  10  9  II  QMODl 
Q3  7  P  3  QMOD2 
•••SECOND**** 

C7  7  14  438P 
C8  12  13  2120P 
C9  14  15  1062P 
CIO  14  0  624P 
R12  7  12  IOOK 
R13  12  15  100K 

R14  13  14  50K 
R 15  14  0  50K 
RI6  18  40  10K 
R17  13  18  51 
RI8  3  16  10K 
Q4  16  15  18  QMODl 
Q5  13  16  3  QMOD2 
•••THIRD**** 

CU  13  21  638P 
C12  19  20  2020P 
C13  21  22  1010P 
C14  21  0  372P 
R 19  13  19  100K 


TAN  it  al AD1C-2.C 
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R20  19  22  100K 
R21  20  21  50K 
R22  2!  0  100K 
R23  24  40  10K 
R24  20  24  50 
R25  3  23  10K 
Q6  23  22  24  QMOD1 
Q7  20  23  3  QMOD2 
•••FOURTH**** 

CIS  20  25  700P 
C16  25  26  700P 
C17  26  29  4000P 
C18  27  28  60000P 
C19  28  0  5000P 
R26  25  0  100K 
R27  25  28  100K 
R28  26  27  100K 
R29  28  29  91K 
R30  28  0  I  OK 
R31  31  40  10K 
R32  27  31  50 
R33  3  30  10K 
Q8  30  29  31  QM0D1 
Q9  27  30  3  QMOD2 
•••FIFTH**** 

C20  27  32  4000 P 
C21  32  35  4000P 
C22  33  34  6600P 
C23  34  0  1300P 
R34  27  34  100K 
R35  32  33  18K 
R36  34  35  10K 
R37  34  0  50K 
R38  37  40  10K 
R39  33  37  50 
R40  3  36  10K 
Q10  36  35  37  QMOD1 
Qll  33  36  3  QMOD2 
C24  33  38  30U 
RL  38  0  200K 


.MODEL  QMOD1  NPN  BF  =  80  VAF  =  200 
+  CJE  =  IP  CJC  =  1P  TF=0.  IN  TR=  100N 
.MODEL  QMOD2  PNP  BF  =  80  VAF  =  200 
+  CJE  =  1P  CJC=1P  TF=0.1N  TR  =  ION 
V1N  I  0  AC 
VCC  3  0  DC  12 
VEE  40  0  DC  -12 
.WIDTH  IN  =80  OUT  =  80 

■  AC  DEC  10  100  20K 
.PLOT  AC  VDB(38) 

■  MIN  AC  OBJT(VDB(7),  VDB(13). 

+  VDB(20).VDB(27),VDB(38)) 

+W=1  -.52  -.52  -.52  -.52  -.52  -.52 

+  -  .5"’  -.52  -.52  -.52  -.52  -.52 

+  -  .20  -.4.0  -8.3  -14.5  -22.4  —33.1 

+  —40.0  -40.0  -40.0  -40.0  -40.0  -40.0 
+  W  =  1  -.48  -.48  -.48  -.48  -.48  -.48 

+  -.48  -.48  -.48  -.48  -.53  -.70  -.90 

+  -1.2  -1.5  -8.7  -26.0  -47.0  -47.0 
+  -47.0  -47.0  -47  0  -47.0  -47.0  -47.0 
+  W  =  1  -.47  -.47  -.47  -.47  -.47  -.47 

+  -.47  -.47  -.47  -.47  -.50  -.60  -.90 

+  - 1.7  -2.4  -4.0  -50.0  -50.0  -50.0 
+  -50.0  -50.0  +  -50.0  -50.0  -50.0  -50.0 
+  W=  1  -41.3  -41.3  -41.3  —41.3  -41.3 

+  —41.3  -41.3  -16  -10  -6  -3  -2  -1.3 

+  -.77  -.38  -2  4  -52  -52  -52  -52 

+  -52  -52  -52  -52  -52 

+  W  =  t  -52  -52  -52  -52  -52  -52  -52 

+  -21  -4  -.2  -.2  -.2  -.2  -.2  -4  -21 
+  -52  -52  -52  -52  -52  -52  -52  -52 
+  -52  XI=0  1  LT  =  1 30  ORDER  =  3 
.PARA  R4  R5  R6  R7  R8  R12  R13  R19 
+  R14  (10K,  1000k)  R15  (I0K,  2000K)  R20 
+  R21  R22  R26  R27  R28  R30  R34  R37 
+  R29  (10K,  2500K)  R35  R36(3K.  2000K) 

■  END 


(5)  An  Inverter  Chain 

CIRCUIT  OF  INVERTER  CHAIN 

.WIDTH  IN  =  80  OUT  =  80 

Ml  3  2  1  1  MODI  L  =20U  W  =  6U 

+  AS=40P  AD  =  20P 

M2  3  2  0  0  MOD2  L  =  6U  W=16U 

+  AS  =  30P  AD  =  30P 

.MODEL  MOD2  NMOS(LEVEL  =  2  VTO=  +  1.5 
+  KP  =  0.8E-5  LAMBDA  =  0.02  CGSO  =  2E-12 
+  CGDO=2E- 12  CGBO  =  3E-  1 1  NSUB=lEi5 
+  C'i  =  2E-7  TOX  =  8E-8  XJ=1U  LD^0.5U) 

CL1  3  0  5E-14 
CL2  4  0  5E-14 

M3  4  3  1  1  MODI  L  =  20U  W  =  6U 

+  AS=40P  AD  =  20P 

M4  4  3  0  0  MOD2  L  =  6U  W  =  16U 

+  AS  =  30P  AD  =  30P 

M5  5  4  1  1  MODI  L=20U  W  =  6U 

+  AS  =  40P  AD  =  20P 

M6  5  4  0  0  MOD2  L  =  6U  W  =  16U 

+  AS  =  30P  AD  =  30P 

M7  6  5  1  1  MODI  L=20U  W  =  6U 

+  AS=40P  AD  =  20P 

M8  6  5  0  0  MOD2  L  =  6U  W  =  16U 

+  AS  =  30P  AD  =  30P 

M9  7  6  1  1  MODI  L=20U  W=6U 

+  AS  =  40P  AD  =  20P 

M10  7  6  0  0  MOD2  L  =  6U  W  =  16U 

+  AS  =  30P  AD  =  30P 

.MODEL  MODI  PMOS  (LEBEL  =  2  VTO  =  -1.5 
+  KP=0.4E  — 5  LAMBDA =0.02  CGSO  =  2E-12 
+  CGDO  =  2E- 12  CGBO  =  3E-  1 1  NSUB=1EI5 
+  a=2E-7  TOX  =  8E-8  XJ  =  1U  LD=0.5U) 

CL3  5  0  5E-14 
CL4  6  0  5E-14 
CL5  7  0  5E-14 
VDD  1  0  DC  5 

VIN  2  0  PWL(0  0  ION  0  20N  5  200N  5) 

.PLOT  TRAN  V(3)  V(4)  V(5)  V(6)  V(7) 

.TRAN  IONS  200NS 

.MIN  TRAN  OBJT(V(3),  V(4),  V(5),  V(6), 

+  V(7))  XI  =0.3  P=4 
+  W  =  5  5  5  0  0  0  (21  points)  0  0  0  0  0 

+  W  =  4  0  0  5  5  5  (21  points)  5  5  5  5  5 

+  W  =3  5  5  0  0  0  (21  points)  0  0  0  0  0 

+  W  =  2  0  0  5  5  5  (21  points)  5  5  5  5  5 

+  W=1  5500000000000000 
tOOOOO  CONTfPV  (2MW)) 

.PARA  Ml,L(5U,40i  >  M2,W(5U,40U) 

+  M3,L(5U,4PU;  M4,W(5U,40U) 

+  M5,L(5U,  r0U)  M6,W(5U,40U) 

+  M7,L.(5U,40U)  M8,W(5U.40U) 

+  M9,L(5U.40U)  M10,W(5U,40U) 

.END 

Appendix  3 

The  Extension  of  POSM  to  Constrained 
Optimization  Applications 

The  optimization  problems  encountered  in  the  field  of 
circuit  and  system  design  are  generally  cataloged  as  the 
constrained  optimization  type.  By  using  the  Lagrange 
multiplier  method  (the  PHR  algorithm  in  [10]),  a  con¬ 
strained  optimization  problem  can  always  be  decomposed 
to  a  series  of  unconstrained  problems. 

The  following  algorithm  summarizes  the  proposed  pro¬ 
cedures  to  deal  with  the  constrained  optimization  prob¬ 
lems  by  employing  the  PHR  method. 

(1)  Select  a  nondecreased  sequence  of  {/***'}  of  pen¬ 
alty  function  factors  and  initialize  tiie  process  by  setting 
initial  point  x°,  initial  multiplier  X(l),  error  tolerance  e, 
and  it  =  1 . 
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(2)  Use  the  proposed  unconstrained  optimization  al¬ 
gorithm  POSM2  to  solve  the  decomposed  unconstrained 
optimization  problem  ((A2)  in  Appendix  1). 

(3)  Update  multiplier  by 

Xr,)-C  +  M,*,«na *{*,■(**). 

=  max  {0,  X}*’  + 

i  =  1.2.  •  •  •  ,/.  (A3) 

(4)  IF 

£  [max  {L,(xk),  <  e2.  (A4) 

THEN  STOP  ( xk  is  the  approximation  solution  of 
the  problem ) 

ELSE  k  =  k  +  1,  GOTO  (2). 


Appendix  4 

The  Derivation  of  S(x) 

Let  Ax,  —  x,  -  xf1,  and  based  on  the  second-order 
error  function  approximation  the  objective  function  can 
be  written  in  implicit  sum  notation  as 

S(X)  =  £  («;(x)]' 

J-  I 
m 

~  £  {ef  +  bn  Ax,  +  Cy Ax2)P. 
jm  1 


This  equation  can  be  expanded  into  binomial  series  like 
S(x)  =  .£  +  p(ef)p',(h,i Ax,  +  c(j4x?) 


,  p(p  0  /  .M\P~2 
+  - jj - (*j  ) 


S(x)  =  .£  (eyM  +  AXj  +  ctj Axj)p 


s  Jw  +  — 
3x,- 


i  a2s 

Ax,-  +  r  t — - — 
ax-o  1  3x,ax,- 


Ax,Ax,- 


3S  3e, 

=  5  + - - 

3e;  3x,- 

32S  3e;  3e; 
3c2  3x,-  3x,- 


Ax  **  0 

,2  . 


Ax  +  if—  d e> 

ax-o  +  2  \3e;  3x,-3x,- 


Ax, -Ax,- 


Ax  “0 


=  5M  +  .£  pePj  \bn  +  2coAx,-)|Ax,0Ax,- 

+  5y?i  M",^y«n-L.0^rAx<. 

1  m 

+  2j?ip(p  ~  eJ  ~2(br,  +  2c,vAx,-) 

•  {bij  +  2coAx,-)|Ax_0Ax,.Ax,. 


=  SM  +  pAxrBR  +  0.5 p{p  -  1) 

•  AxTBDBTAx  +  PYTCR  ( A6) 

where  5,-,-  is  Dirac  5  function.  (A5)  is  identical  to  (A6) 
except  for  the  inclusion  of  two  additional  high  order  terms 
which  usually  do  not  introduce  too  much  difference  and 
little  improvement  on  approximation.  Thus  we  adopted 
(A5)  rather  than  (A6)  in  this  algorithm  for  more  accurate 
approximation. 
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■  {bh  Ax,  +  c,j\x?)2  +  •  •  •  j 

-s"+£(  {,(*-)'-■  ■■f.  +  pv^rVx? 

+  0.5 p(p  -  1  )(e  'byAXibfjAx,' 

+  0.5 p(p  -  l)(e"/  xfcijAxf- 

+  P(P  ~  !)(<V  ~b,jAx,crjAxl  +  •  •  •  }• 

Neglecting  the  higher  order  terms,  we  can  express  S(x) 
in  a  more  compact  form  as 

5(x)  ~  Sw  +  pAxTBR  +  0.5 p{p  -  \)AxTBDBT Ax 

+  pYTCR  +  0.5 p(p  -  1  )YtCDCtY 

+  p{p  -  \)AxtBDCtY.  ( A5) 

If  the  Taylor  series  expansion  is  assumed,  the  objective 
function  will  be 
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